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ABSTRACT 

The performance of an unaided attitude determination system based on GPS interfer- 
ometry is examined using linear covariance analysis. The modelled system includes 
four GPS antennae onboard a gravity gradient stabilized spacecraft, specifically the 
Air Force’s RADCAL satellite. The principal error sources are identified and mod- 
elled. The optimal system’s sensitivities to these error sources are examined through 
an error budget and by varying system parameters. The effects of two satellite selec- 
tion algorithms, Geometric and Attitude Dilution Of Precision (GDOP and ADOP, 
respectively) are examined. The attitude performance of two optimal-suboptimal fil- 
ters is also presented. Based on this analysis, the limiting factors in attitude accuracy 
are the knowledge of the relative antenna locations, the electrical path lengths from the 
antennae to the receiver, and the multipath environment. The performance of the sys- 
tem is found to be fairly insensitive to torque errors, orbital inclination, and the two 
satellite geometry figures-of-merit tested. 
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Chapter 1 


Introduction 

This introductory chapter first provides the motivation for and a brief background 
on the use of the Global Positioning System (GPS) in attitude determination. The 
objectives of this investigation are then outlined, followed by an introduction to the 
basics of GPS interferometry and the types of measurements which can be made using 
the interferometer to solve for three-axis attitude. Then follows a description of the 
interferometer considered in this investigation: four patch antennae onboard a gravity 
gradient stabilized spacecraft modelled after the RADCAL satellite. 

1.1 Background 

GPS technology provides a solid-state continual solution to the problem of solving 
for attitude, previously restricted to gyro-based systems. This technology offers high 
reliability, very little maintenance, and lower costs while potentially achieving millira- 
dian-level accuracy [21]. Attitude determination using GPS is based on a simple 
geometric principle: two distinct points uniquely define a line and three noncolinear 
points define a plane. Treating GPS antennae, with known relative positions, as the 
noncolinear points, the orientation of the plane in which they reside can be deter- 
mined, as well as the orientation of the antennae within the plane. This is achieved by 
using the differences in phase of the incoming GPS carrier signal between the anten- 
nae, hence the term “GPS interferometry.” 

The use of GPS carrier phase difference measurements for attitude determination 
of a stationary platform was first reported by Texas Instruments in 1981 [21]. Other 
manufacturers including Trimble Navigation, Ashtech, Adroit Systems, and Loral 
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have also developed GPS receivers for attitude determination. The capability of GPS- 
based attitude determination has been demonstrated in land-, air-, and marine-based 
applications [6, 19, 21, 27]. The first flight tests of a real-time GPS attitude determina- 
tion system for an aircraft were conducted at Ohio University in 1991 [27]. Currently, 
there is no known use of GPS interferometry in a closed loop system. 

The application of this method of attitude determination has recently been 
extended to space-based platforms. The GPS Attitude and Navigation Experiment 
(GANE), which will test a GPS interferometer intended for International Space Station 
Alpha (ISSA), is currently scheduled to fly onboard the Shuttle in April 1996 [9]. 
Mounted on a 1.5 by 3 m platform in the cargo bay will be a four-antenna interferome- 
ter, along with an Inertial Reference Unit (IRU) for attitude verification. The require- 
ments for this stand-alone GPS interferometer are to estimate the station’s attitude to 
within 0.3 deg (3a) per axis and attitude rates to within 0.01 deg/s (3a) per axis at 0.5 
Hz. 

An ongoing spacebome experiment in GPS interferometric attitude determination 
is flying on the Air Force Space Test Program’s RADar CALibration (RADCAL) sat- 
ellite. RADCAL is the specific application of interest in this investigation. The RAD- 
CAL satellite was launched from Vandenberg Air Force Base on June 25, 1993. Its 
principal purpose is for calibration of civilian and military C-band radar stations. The 
spacecraft payload includes a radar transponder, a TRANET beacon for navigation, a 
magnetometer, and a Trimble Advanced Navigation Sensor (TANS) Quadrex receiver, 
specially adapted for attitude determination by Stanford University [8]. The space- 
craft is gravity gradient stabilized, but no other source of attitude information exists on 
RADCAL, save the magnetometer. 
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Several research groups are involved in the RADCAL experiment in attitude 
determination, including Cohen, Lightsey, and Parkinson at Stanford University, Axel- 
rad and Ward at the University of Colorado, and others at JSC. Cohen, et al, report a 
post-processed attitude solution accuracy (la) on the order of one degree, though they 
ultimately expect to achieve an accuracy of 0.3 deg in the post-flight solution [8]. 
Axelrad and Ward are working to implement their “bootstrapping” algorithms, 
intended for onboard navigation, to provide high accuracy (~ 0.2 deg/axis) attitude 
solutions for near Earth satellites [2, 28] and are testing these algorithms on RADCAL 
flight data. The bootstrap process starts with an initial attitude estimator which 
resolves the integer ambiguities and approximates line biases, followed by a baseline 
estimator, and finally a Kalman filter which provides the best attitude estimate in real- 
time. 


1.2 Objectives of This Investigation 

Developed in the following pages is a linear covariance analysis of a GPS interfer- 
ometric attitude determination system for a gravity gradient stabilized spacecraft, 
using the RADCAL satellite configuration as an example. The goals of this study are 
to determine how well attitude can be maintained by processing GPS measurements 
with a Kalman filter and to understand the nature of the error sources which compro- 
mise the attitude solution. System parameters including noise levels, satellite selection 
criterion, and orbital inclination are varied to gain insight into the optimal filter’s per- 
formance. An error budget is tabulated and used in determining the major contributors 
to the attitude rms errors. In anticipation of designing an onboard filter, the perfor- 
mance of a reduced-order or suboptimal filter is analyzed using consider states. Unfor- 
tunately, it was not possible to process actual RADCAL flight data in the form of raw 
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single phase differences, obtained from Johnson Space Center, in the time available. 
Contained in an appendix is a discussion of this data and some of the problems associ- 
ated with it. 

1.3 GPS Interferometry and Observables 

The primary observable in attitude determination using GPS is the fractional dif- 
ference in carrier phase, A<j), between two antennae. The carrier signal is normally at 
•the LI frequency for maximum accuracy, centered on 1575.42 MHz, corresponding to 
a wavelength of roughly 0. 19 m. For multiple antennae, the phase differences are typi- 
cally referred to a common antenna, identified as the master, and the others are desig- 
nated as slaves. Figure 1.1 depicts the differential phase geometry of a GPS 
interferometer consisting of a single master-slave antenna pair. The phase difference is 
related to the difference in range, A r (in cycles), from the satellite to each antenna via: 

A(j) = Ar-n + p + £ + v (1.1) 

where n is the integer number of cycles in the differential range, f3 is the line bias, £, the 
multipath, and v the receiver noise. Several techniques with varying degrees of com- 
plexity exist for determining the integer ambiguity, but it need only be resolved once, 
barring future cycle slips. The line bias given here is actually the difference in path 
delays from each antenna to the receiver. As can be seen from Figure 1.1, the range 
difference is merely the component of the baseline, b, projected onto the unit line-of- 
sight (LOS) vector, e, and is related to the attitude as follows: 

A r = b r e = £cos0 (1.2) 
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Figure 1.1 GPS Interferometer Geometry 


Resolving the baseline and LOS vectors into a common frame is equivalent to deter- 
mining the attitude of the baseline: 

Ar = b g M 5£ e£ ( 1 . 3 ) 

where M B £ is the transformation from the frame in which the line-of-sight vector is 
expressed, usually an Earth-Centered-Earth-Fixed (ECEF) frame, to that of the base- 
line vector, usually a body-fixed frame. 

Note that any rotation about the baseline axis is not observable, thus requiring a 
second baseline to determine the host body’s three-axis attitude. Note further that a 
rotation about the LOS vector is also unobservable. Therefore, any number of base- 
lines with just one satellite is insufficient to determine M BE . Expressing position only 
requires three pieces of information. For full attitude specification, however, this is not 
quite adequate. In general, attitude specification requires three Euler angles and the 


sequence. This fact holds true even if the attitude is expressed as a quaternion or trans- 
formation matrix. However, only small corrections about a nominal attitude are dealt 
with here, and these corrections can be applied in any order. Therefore, only three 
pieces of information, i.e., measurements, are necessary to determine attitude errors 
(though two baselines are still required). 

An attitude correction requires one less measurement than the four required for 
translational error state (three position, one clock) determination. This is a conse- 
quence of the differential nature of interferometry. The differencing operation is done 
with respect to a common oscillator in the receiver hardware, thus effectively eliminat- 
ing user clock error. The single difference is free of the offset and the linear drift of the 
satellite frequency [17]. Small errors due to the effect of the changing satellite fre- 
quency between signal transmission times for the two antennae are present, but due to 
the proximity of the antennae and the high quality of the satellite oscillator, these 
effects are negligible. The ionospheric delay experienced by each signal is also essen- 
tially the same and cancels out in the single difference. These and other errors which 
corrupt the attitude solution will be discussed later in more detail. 

Though only the single phase difference measurement is used in this investigation, 
two other combinations of the carrier phase measurement can be used in interferomet- 
ric attitude determination: the double difference and the triple difference [17, 29]. The 
double difference is typically formed by differencing phase difference measurements 
from two satellites, though the difference can be a temporal one. In the between-satel- 
lite double difference, the error resulting from the difference in electrical path length 
to each antenna is removed. The double difference has its disadvantages, though, since 
the receiver noise is magnified in combining the phase difference measurements, and 
another satellite is required to obtain the same number of measurements. This might 
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pose a problem if the antennae do not view the same region of sky. The triple differ- 
ence is formed by differencing over antennae, satellites, and time. The advantage of 
this observable is that the initial integer ambiguity, which is constant in time, is elimi- 
nated (assuming no cycle slips between epochs). 

1.4 Interferometer Configuration 

The interferometer in this investigation resides on a vehicle modelled after the 
RADCAL satellite, which is in a polar, near-circular, 815 km altitude orbit. It is grav- 
ity gradient stabilized, requiring the incorporation of a dynamic model in the attitude 
determination process. Figure 1.2 depicts the RADCAL satellite and the GPS antenna 
array. Four patch antennae are mounted on the zenith face of the axially symmetric 




Figure 1.2 RADCAL Configuration 

craft and are equally spaced about the perimeter of the 30-inch diameter cylindrical 
bus. Each antenna is canted away from the zenith direction by 17.5 deg in order to 
maximize GPS satellite visibility while reducing multipath from the magnetometer 
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boom. The receiver time-multiplexes through the four antennae. It simultaneously 
closes 24 separate tracking loops and can therefore provide phase measurements to 
each of the four antennae for up to six satellites per antenna, or up to 18 phase differ- 
ence measurements per epoch for the three baselines [8]. 

The Quadrex receiver onboard RADCAL is overly restrictive in that it records 
phase difference measurements only to satellites in view of all four antennae, even 
though a measurement could be processed if only two antennae view the satellite. The 
system is nevertheless still overdetermined given the 18 measurements to the three 
baselines. The simulation developed here employs the phase differences over three 
baselines with antenna 1 as the master. The biases from the oscillator and receiver 
electronics drop out in forming the phase differences, resulting in measurement noise 
errors of less than a millimeter [8], A greater contributor to the error is multipath, the 
effect of which is speculated to be about 5 mm for the RADCAL configuration [8, 18]. 
These and other error sources which corrupt the attitude solution are described further 
in the following chapter. 
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Chapter 2 


Error Modelling 

A survey of the literature reveals that there are many error sources which can limit 
the attitude solution accuracy using GPS interferometry. This chapter outlines a few of 
the major contributors and discusses the modeling of those used in this analysis. 

2.1 Receiver Noise 

The calculation of the incoming signal’s phase angle is subject to error induced by 
the hardware itself. The phase error depends on the receiver’s quality, as well as the 
vehicle dynamics. The dynamics of the gravity gradient stabilized satellite are slow 
enough such that their risk of not being within the tracking loop bandwidth is negligi- 
ble. The Quadrex receiver has 24 tracking loops and is capable of an output rate of 10 
Hz, for a total of 240 phase measurements per second. According to Cohen, Lightsey, 
and Parkinson, the Quadrex’s architecture provides very clean phase measurements 
[8]. As previously mentioned, the biases from the frequency standard and receiver 
hardware cancel in differencing these measurements. Thus, phase difference measure- 
ments should be corrupted by less than a millimeter after passing through the RF 
switch at the receiver’s front-end. 

2.2 Vehicle Flexibility 

Knowledge of the spacecraft attitude can only be as good as that of the baseline 
vectors, for it is these vectors which define the orientation. A trade-off exists in choos- 
ing the baseline length. One might expect that greater angular resolution would be 
obtained with longer baselines, but more integer ambiguity and unaccounted vehicle 
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flexing between the antennae might negate that benefit [15]. An acknowledgment that 
the baselines are subject to flexures is therefore required for realistic attitude determi- 
nation. One source of flexure would be thermal stresses experienced by the craft when 
passing in and out of eclipse. The effect of the flexure can be to change the baseline 
length, as well as its direction. The latter result may not be easily differentiated from a 
change in vehicle attitude. Flexures of the RADCAL baselines are kept reasonably 
small since the antennae are placed on the fairly rigid platform of the spacecraft bus 
and are less than a meter apart. 

2.3 Line Bias 

The line bias over a baseline results from the difference in electrical path lengths 
from each antenna to the receiver. This “bias” may not be constant, however, since the 
effective path lengths may change with thermal variations, similar to the baseline vari- 
ations previously mentioned. Minimizing the error can be accomplished by keeping 
the path lengths as short and as symmetric as possible, and by configuring the vehicle 
such that the cabling is subject to small temperature gradients [3, 6]. Fortunately, the 
phenomenon affecting the path delay to one antenna should be similar to that affecting 
the path delay to another; the changes in path delay are therefore mitigated somewhat. 
The line bias can be eliminated by utilizing the between-satellite double difference 
observable. In the sequel, the terms “path delay” and “line bias” will be used inter- 
changeably to denote the path delay to one antenna, unless specified otherwise. 

2.4 Multipath 

Multipath is the undesired reflection of the incoming GPS signal from the 
antenna’s surroundings. The antenna may receive both the direct and reflected signal, 
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or could conceivably receive only the reflected one. Since the reflected signal travels a 
different path, its phase is shifted from the direct signal. The reflected signal then 
appears as an additive bias to the primary transmission, as indicated in Equation (1.1). 
For small baselines on a spacecraft, the multipath error experienced by each antenna 
could have a common component. In the RADCAL scenario, however, the antennae 
are canted outward and surround the base of a 6-meter boom; the majority of the mul- 
tipath is then likely differential mode and does not drop out in the phase difference. 

Multipath error can be diminished by various techniques. A low-multipath envi- 
ronment is most desirable; antennae, therefore, should obviously not be placed near 
multipath sources, though this is unavoidable in the RADCAL configuration. A suit- 
able choice of coating on the mounting surface can also reduce reflectivity [21]. The 
antenna gain pattern can be appropriately shaped to mask out signals entering from 
directions of suspected multipath sources. Another way to mitigate multipath, 
described in great detail by Cohen and Parkinson, is calibrating out the repeatable part 
of the effect [6, 7]. The only source of multipath in a spacecraft application is the 
spacecraft itself. A wave front from a given direction will reflect off vehicle surfaces 
in a repeatable way, provided that the reflectivity does not change and there are no 
moving parts in the viewing environment. The multipath can then, in theory, be cali- 
brated out as function of incidence direction, and whatever error remains is receiver 
noise. Cohen goes on to develop a spherical harmonic model of the repeatable multi- 
path characteristic for an experimental configuration. He suggests that the residual car- 
rier phase error can be reduced from 5 mm for the RADCAL configuration to less than 
1 mm with multipath calibration [8]. 
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2.5 Antenna Phase Center Variations 


Another error source is the asymmetry of the antenna gain pattern [1, 18]. This 
asymmetry can cause significant differences in phase based on the signal’s incidence 
direction, resulting in an apparent migration of the phase center. Two antennae of the 
same make can have similar gain patterns, and the error can be small if the two are 
similarly oriented with respect to the incoming signal. But if the signal arrives from 
two different antenna-relative directions caused by, say, the antennae being canted 
away from each other, the error in differencing the two phases may not be so small. 
Cohen asserts that this effect can be combined with multipath if one considers the 
entire vehicle to be the antenna, and due to their dependence on the direction of signal 
incidence only, the two error sources can be calibrated out [6, 7]. 

2.6 Dilution of Precision 

The attitude solution is also compromised by the GPS satellite geometry, itself 
[18]. For translational space solutions, the familiar Geometric Dilution Of Precision 
(GDOP) figure-of-merit is used to assess the impact of satellite geometry on position 
and time determination. A smaller GDOP corresponds to a greater volume of the poly- 
hedron formed by connecting the vertices of the unit vectors from the user to each 
GPS satellite. But for attitude determination, the best resolution occurs when the line 
of sight is perpendicular to the baseline vector; an alternate figure of merit, Attitude 
Dilution Of Precision (ADOP) [26], is therefore required to rate the effect of satellite 
geometry on the attitude solution. (The derivations of these criteria can be found in 
Section 4.2.) Note that with only one satellite visible, rotations about the LOS are 
unobservable. In RADCAL’s viewing environment, this situation is not much of a 
problem, though satellite geometry will still affect the attitude solution. 
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2.7 Integer Ambiguity and Cycle Slip 


Obviously, the GPS receiver cannot distinguish one cycle of the carrier signal from 
any other. Thus, for baselines longer than half of a wavelength (LI carrier wavelength 
= 19 cm), an ambiguity in the integer number of cycles between the two antennae 
exists [15]. Translating the phase difference, A<f>, into the range difference, A r, and ulti- 
mately solving for the relative positions of the antennae requires resolution of this 
integer ambiguity. The concept of the ambiguity is illustrated in the two-dimensional 
example in Figure 2.1, for one satellite and aim baseline [15], The master antenna is 
at the center, and the slave resides somewhere on the 1 m radius circle. The vertical 
lines are the lines of constant phase. The range difference contains the unknown inte- 
ger, n, and the measured (noiseless) phase difference for this example is A<J) = 0.0. As 
one can see, there are 1 1 possible range differences (and n’s) to give 22 possible base- 
line orientations, i.e., attitudes. 



Figure 2.1 Integer Ambiguity for lm Baseline & One Satellite 


In the three-dimensional problem, still with one baseline and one satellite, the set 
of possible attitudes defines a sphere with a 1 m radius, and the ambiguities define a 


set of planes which are perpendicular to the LOS and intersect the sphere. The loci of 
intersections, in other words, are a family of parallel circles whose center is on the 
LOS. A second satellite creates another family of circles. The number of possible solu- 
tions is now reduced to the points of intersection between the two families of circles. 
In reality, however, the measurements are noisy - the points of intersection then 
become regions of uncertainty. The addition of more satellites will reduce the number 
of possible solutions. The added constraint introduced with another antenna can also 
be used to limit the search, by utilizing the baseline relative geometry. 

For the three baselines which view up to six satellites at a time (typical of RAD- 
CAL), there will be 18 integer ambiguities. Once the ambiguities are resolved, they 
can be maintained until a cycle slip occurs, resulting from a loss of lock by the track- 
ing loop. If the attitude prior to the slip is known well enough, the slip can be detected 
and resolved. For the longest RADCAL baseline, 62.6 cm, attitude accuracy better 
than 8 degrees is, in theory, sufficient to uniquely determine the integer ambiguity. 

Many algorithms to resolve the integer ambiguity have been suggested with vary- 
ing degrees of efficacy. One technique is the brute force integer search, in which dif- 
ferent combinations of integers are checked. Methodical searches which utilize 
geometric constraints to limit the search have been proposed [4, 11, 14]. Knight 
describes additional time saving techniques in “pruning the decision tree” in the search 
and claims to resolve the integers instantaneously (using measurements at one epoch) 
for a four antenna array on a two-meter square with four satellites in track with 99.3% 
success rate in an average time of 0.24 seconds on a 486 PC [14], The success rate 
may not be so great in a high-multipath environment. The problem with the integer 
search is that it may conclude with an incorrect result, since there may not be a unique 
solution [6], depending on the baseline geometry and the noise level. An alternative to 
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the integer search suggested by Cohen is a motion-based approach, which makes use 
of the additional information provided by the baseline motion [6]. The process 
involves mapping the phase differences into a set of baseline displacement vectors. 
Assuming the baseline length to be constant, and no cycle slips from one measurement 
to the next, the baseline solution may be obtained through a linear least-squares fit, for 
large angle motion. This solution is used as an initial guess and is refined using a non- 
linear least squares fit to best match the raw phase measurements. With an initial atti- 
tude, the integer ambiguities fall out. Cohen also outlines a scheme for small angle 
motion. Whatever the method, the ambiguity must be reconciled to determine attitude 
accurately. 

2.8 GPS SV and User Clocks 

The clocks of both the GPS satellite and the user receiver are always in error to 
some extent. In forming the phase difference observable, however, the majority of that 
error is removed. The large terms resulting from the offset and linear drift of the GPS 
satellite frequency cancel [17]. The smaller terms from the errors in satellite frequency 
between the signal transmission times are negligible due to the proximity of the anten- 
nae and the stability of the GPS clocks. If the receiver position is determined from the 
GPS Coarse/ Acquisition (C/A) code, the solution will be subject to the satellite clock 
dithering effect known as Selective Availability (SA). Virtually no accuracy is lost in 
the attitude solution due to SA, however, since the lines-of-sight are insensitive to the 
position errors induced by the intentional degradation. The user clock errors cancel 
altogether, since the signals from all antennae are referred to a common oscillator. 
Thus, clock errors tend not to affect attitude determination. 
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2.9 Propagation Media 

The propagation media affect electromagnetic waves at all frequencies, resulting 
in refraction of the rays which manifests itself as a time delay of arriving signals. The 
troposphere, up to 40 km altitude, is a nondispersive medium, in which the refraction 
is independent of the GPS signal frequency, and it imposes an apparent increase in 
path length [17]. The ionosphere, from 100 to 1000 km, in which RADCAL resides, 
has an effect which depends on the total electron current (TEC) or the total number of 
electrons in the path between GPS satellite and receiver. This effect, arising from the 
dispersive nature of the medium, is to retard the modulation on the carrier wave while 
advancing the carrier phase [16]. Fortunately, the signals arriving at the interferometer 
travel nearly identical paths, and any delay or advance is effectively eliminated in dif- 
ferencing phase measurements. (The cabling from each antennae to the receiver is also 
dispersive, indicating that the path delays mentioned earlier could actually advance the 
carrier phase [12]. This dispersion, unlike that of the ionosphere, is likely to be nearly 
constant and is included in the line bias.) 

2.10 Modelling of Error Sources 

Several contributors from among the previous list of error sources are character- 
ized in this investigation using statistical models. These are: receiver noise, baseline 
motion (flexibility), line bias, and multipath. The receiver noise, or total phase mea- 
surement error for each antenna is modelled as zero-mean, Gaussian white noise. 
Phase errors between antennae, and between satellites are assumed to be initially 
uncorrelated. Correlations do arise, however, when the phase difference measurements 
are formed. This topic will be discussed in more detail in Chapter 4. 
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The effects of baseline flexures and line biases are included in this analysis even 
though they are not explicitly modelled. Rather, each antenna’s lever arm (from the 
center of mass) flexure and electronic path delay error state are modelled as first-order 
Markov processes. The flexure or delay becomes exponentially decorrelated over time 
as described by the autocorrelation function, R x , of the Markov process [5, 10]: 

H M 

R x (At) = a 2 e T ( 2 . 1 ) 

where a 2 is the mean-square value of the process, and t is its time constant, or the time 
separation required for the correlation to be reduced by 63%. In choosing an appropri- 
ate time constant, we must consider the physical phenomenon at work. In this case, we 
assume the driving force of the errors is thermal variation caused by solar heating on 
the rotating spacecraft. It therefore seems reasonable to choose a time constant of half 
an orbital period. A more sophisticated and accurate model of the flexure dynamics 
could be used, but the Markov description should adequately characterize the process. 

Multipath is also modelled. Without a spherical harmonic model of the multipath 
environment for RADCAL, however, the simpler Markov model is assumed. Cohen 
[8] suggests that a mean-square value of 5 mm is appropriate. The time constant is the 
time required for the LOS to the satellite in view to ‘significantly’ change direction, 
defined as a rotation of 30 degrees [19, 26] in the body frame. The multipath is 
assumed to be independent from one satellite to another, though from any one satellite, 
the multipath may be correlated among the antennae. The degree of correlation 
depends on the environment’s (i.e., vehicle’s) local reflectivity and is quantified by a 
variable correlation coefficient, p^. See Chapter 4 for more detail. The actual multipath 
behavior is more complicated, but this model can be used as a simple means to 
observe its effects. 
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As indicated earlier, clock and propagation media errors are considered negligible 
because of the differential nature of the problem. The antenna phase center variations 
can be lumped in with the multipath, since the two effects are practically indistinguish- 
able. The error arising from poor satellite geometry is tempered by employing a satel- 
lite selection algorithm. The two algorithms considered are based on GDOP and 
ADOP. GDOP is used since RADCAL selects satellites based on this criterion for its 
translational state estimation. ADOP is used for comparison purposes because it more 
correctly assesses the impact of satellite geometry on the attitude solution. Note that 
the two algorithms may cause different groups of satellites to be selected. 

The integer ambiguity is assumed known in the linear covariance analysis. Of 
course, in filtering real flight data, the ambiguity would have to be resolved at least in 
determining the initial attitude. Measurements were assumed to be cycle-slip free for 
the covariance analysis, because, although real data will contain cycle slips, they are 
assumed detectable. 
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Chapter 3 


Attitude Error Dynamics 

The reference frames used in this analysis are described herein, followed by an 
analytic development of the attitude error dynamics for gravity gradient stabilized 
spacecraft. Also presented is the analytical solution for the specific case in which the 
nominal attitude is coincident with the local-vertical-local-horizontal (LVLH) frame. 

3.1 Reference Frames 

Four reference frames are employed in this investigation. The first is the Earth- 
Centered-Earth-Fixed (E) frame in which the GPS SV positions found in the Interna- 
tional GPS Service (IGS) precise ephemeris files are expressed. (These precise ephe- 
merides are expressed in the International Terrestrial Reference Frame (ITRF-93).)The 
z axis is aligned approximately with the spin axis of the Earth, the x axis lies in the 
equatorial plane and points toward the Greenwich Meridian, and the y axis completes 
the right-handed orthogonal set. An Earth-Centered-Inertial (I) frame is defined as the 
E frame at time zero - the transformation from the I to E frame is merely a rotation 
about the inertial z axis (see Figure 3.1). This approximation is valid only for short 
periods of time (on the order of days) due to the drift of the Earth’s spin axis. 

The third reference frame utilized is the LVLH (L) frame, centered on the vehicle. 
The z L vector points in the negative radial direction, the y L vector points in the nega- 
tive orbital angular velocity direction, and the x L vector is horizontal forward (see Fig- 
ure 3.2). 

Finally, the Body-Fixed (B) frame is aligned with the spacecraft’s principal axes 
and nominally coincides with the L frame due to the gravity gradient stabilization and 
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RADCAL’s axial symmetry (see Figure 3.2). The orientation of the B frame, and thus 
the body attitude, is defined by the successive Euler angle rotations: first yaw (0 3 ) 
about the z L axis, then pitch (0 2 ) about the new y axis, and finally roll (0]) about the 
resulting x axis, now \ B . The transformation from the L frame to the B frame can be 
written: 


= T Jt (0 1 )T y (0 2 )T z (0 3 ) 


where 



1 0 0 


cos0 2 0 -sin0 2 

T ,(6,) = 

0 cos0, sin0, 
0 -sin0, cos©. 

T y (0 2 ) = 

0 1 0 
sin0 2 0 cos0 2 


T z (9 3 ) 


cos0 3 sin0 3 0 
-sin0 3 cos0 3 0 
0 0 1 


The above attitude matrix can be linearized about an estimated (or nominal) atti- 
tude in the following manner: 

M, t -M, 4 M <1 -(I-[Sex])M ta (3.,) 

where is a transformation from the estimated body frame, B , to the true body 
frame, B, through the small angles in the vector, 50: 


50 = 0-0 where 0 = 


roll 

pitch 


yaw 


and the associated cross-product matrix is: 


[ 60 x] 


0 -50 3 50 2 
50 3 0 -50, 

-50 2 60, 0 


Note that rotations through small angles are commutative. 


( 3 . 2 ) 
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Figure 3.1 Inertial, Earth-Fixed, and Local Coordinate Frames 



Figure 3.2 Local and Body Coordinate Frames 

3.2 Attitude Error Dynamics of a Gravity Gradient Stabilized 
Spacecraft 

The focus of the linear covariance analysis is the error statistics’ evolution, which 
depends on the error dynamics. The goal in this section is to derive the dynamics of 
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the errors given in Equation (3.2) and of their time derivatives. The starting point is 
Euler’s vector equation of motion, expressed in a body-fixed frame: 

(CO) = I ^ 1 (T - COX (I M CO) ) (3.3) 

where 

co = angular velocity vector, expressed in body coordinates 
l M = moment of inertia matrix expressed in body coordinates 
T = external torque on the rigid body expressed in body coordinates 

3.2.1 The General Case 

In rewriting Equation (3.3) in terms of the Euler angles and their errors from the 
nominal body frame, the following expressions will be required: 

^(Myg) = M /b [cox] (3.4) 

5(0 = (0- ft) (3 5) 

The first is an expression for the time derivative of the attitude matrix, M /JS , which is 
the transformation matrix from the body frame to the inertial frame. The second is the 
angular velocity vector error. Let us assume that there are no translational errors nor 
inertia matrix errors; that is, errors only arise from the uncertainty in the transforma- 
tion to the true body frame. The derivation of the sensitivities proceeds [24], 

Given the above definitions, we can express the transformation from the true body 
frame to the inertial frame in terms of the transformation from the nominal body frame 
and the attitude error, to first order: 

M„ = M /A M to = M /B (I + [80x] ) (3.6) 

Substitute this expression and Equation (3.5) into Equation (3.4): 
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(3.7) 



or 

4.(89) = 50 = 8co- [d)x] 56 (3.10) 

at 

Take the time derivative of both sides: 

^(59) = ^(8co) - ( co) xj 50- [ tbx] ^ ( 50) (3.11) 

The first term on the right can be rewritten by linearizing Equation (3.3): 

^(®) =^(« + 5to) = I M ‘ [T + 6T- (cb + 8co)x(I M (G) + 5a)))] (3.12) 

yielding the first-order perturbation equation: 

^(5(0) = i;'[§T+ [[(I»x] - [c bx]I„] (8(0)] (3.13) 

The torque error can be decomposed into the contributions from the gravity gradient 
and other sources: 


8T 


^ cc 
39 


80 + 


—^80 + 8T 

39 


other 


(3.14) 


The second term on the right is zero since the gravity gradient torque depends only on 
the attitude and not its rate, as will be seen. The last term in the above equation embod- 


35 


ies all other sources of torques: solar radiation pressure, atmospheric drag, cosmic 
dust, and other perturbations. The modelled spacecraft has no control torques acting 
on it. Now insert equations (3.13), (3.14), and (3.10) into Equation (3.11) and factor 
out 80, Sco, and &T other to get: 

^(50) = l©x] + [fix] [fix]] 86 (3.15) 

+ Pi [ [ (I» x] - [fix] l M ) - [fix]]8co 

+ ^8T other 

Now, substitute 8to from Equation (3.10) and let 8T = ST other : 

4(60) = 80 = A80 + B80 + C8T (3.16) 

at 

where 

r 9 t n 

A = I' 1 [ [ [ (I M fi) x] - [d>x] I M ] [fix] + - [fix] (3.17) 

B = I"‘ [ [ (I» x] - [fix] IJ - [fix] (3.18) 

c = rj (3.i9) 

We still must get an expression for the sensitivity of the gravity gradient torque to the 
attitude in Equation (3.17). The gravity gradient torque is, in body coordinates [13, 
30]: 

T gg = ^[rxtVO] (3.20) 

M 

Since our knowledge of the body frame is in error, the position vector in (nominal) 
body coordinates is in error. Using the correction transformation from the nominal 
body frame to the true body frame given in Equation (3.6), express the true position 
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vector in terms of the attitude error, 50, and the position vector in nominal body frame 
coordinates, r : 

r = M g -f = (1+ [50x] ) r f = (I - [50x])f (3.21) 

Likewise, express the torque as a nominal term plus a correction and substitute the 
above into Equation (3.20) to get: 

T go + 8T co = ^ [fx (I w f) - rx (I M (80xr) ) - (80xr) x (I w f) ] (3.22) 

If | 

again, ignoring terms of second-order and higher. The perturbation in the gravity gra- 
dient torque is therefore: 

5 t g<; = -gf^Se (3.23) 


where 

0T o 

- 5 ^ = it [ [fx] I M - [ (I„f) X] ] [fx] (3.24) 

CTO |fj 

Substitute this result into Equation (3.17), and we now have an expression for the time 
rate of change of the attitude error rate in terms of 50, 50 and 5T. The attitude error 
dynamics equation in matrix form is given as: 


50 

50 


^ 3*3 ^ 3*3 ®3 
1^3*3 ®3jt3 ^3 


-] 

50 


50 

y 

ST 


(3.25) 


where A, B, and C are given in expressions (3.17), (3.18), and (3.19), respectively. 


3.2.2 Special Case 

For the nominal case, in which the small rotations are from the LVLH frame, 
Euler’s equations describing the motion of the gravity gradient stabilized spacecraft 
can be simplified. Let us assume the body frame coincides with the principal axes, the 
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orbit is circular, and only the gravity gradient torque is acting (8T = 0). The inertia 
matrix, nominal angular velocity and position vectors reduce to: 



\ 0 o' 


0 


0 

Im = 

o 

o 

ft) = 

-<»orb 

= const f = 

0 


0 0 I z 


0 


-r 


2 <3 

where (£> orb = fi/r is the orbital rate for a spacecraft in a circular orbit of radius, r, 
about a point mass with a gravitational constant, p. Inserting the above quantities into 
equations (3.17), (3.18), and (3.24), respectively, the linearized equations of motion in 
terms of the differential roll, pitch, and yaw angles are: 


50,+ 4(0^60, -co„ r ^5e 3 = 0 

(3.26) 

502 + 3(0^0^502 = 0 

(3.27) 

503 + (oLo 2 80 3 + co„ r ^5'0i = 0 

(3.28) 


where 




a = -4 — - 2 
x 


h-h 


L-L 




x r ** 7 


The modelled craft is axially symmetric with I* = Iy = 26.4 kg m 2 , and l z = 5.813 kg 
m 2 and resides in an 815 km altitude orbit (the principal moments of inertia and 
approximate altitude of RADCAL [18]). From Equation (3.27), we readily see that the 
pitch is decoupled and is a pure sinusoid of the form: 


59,(0 = n,sin ( ((O orb j3c y ) t+ a 2 ) (3.29) 

where a x is a constant of integration. The pitch is seen to have a period of oscillation of 
2nl ((O orb j3o y ) ~ 3960 s. The differential roll and yaw solutions can be found by first 
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rewriting equations (3.26) and (3.28) in matrix form [18]: 




1 0 0 80, 

0 0 ®orbK 56 , 

0 0 1 80 3 

-®„rb k z 0 803 


Noting that a z = 0 and k, = 1 for this particular case gives the following eigenvalues 
for the roll and yaw modes: 


roll: \ 2 = ±U or J(3o x + 1 )i 

yaw: X 3 4 = 0, 0 


Taking these eigenvalues and their corresponding eigenvectors gives the roll and yaw 


solutions: 


80,(0 = a 3 sin ( (a orb j3a x + l)t + a 4 ) 


80 3 (O = . * 3 cos ( (® ()rb j3o x + 1) f + a 4 ) + a^t + (3.31) 


< 73 °, + 1 

The roll and yaw are seen to have periods of oscillation of 2iz/ (® (>rb j3o x + 1) ~ 
3320 s, though the yaw has a constant offset and a linear term. Equations (3.29) 
through (3.31) may also be found in Reference [18]. This behavior is verified by 
directly integrating equations (3.3) and (3.4). A fourth-order Runge-Kutta scheme was 
used to integrate these two equations for a circular 815 km altitude orbit using a two- 
body gravity model, with the following initial conditions: 


80 (r = 0) = 60, = _ 5 _o deg, to L (t = 0) = -© = -0.0593 deg/s 


39 


where ©l is the orbital angular velocity expressed in the local frame. Figure 3.3 shows 
the resulting roll, pitch, and yaw, which are extracted from M /B , for a four-hour time 
span. Note that for the small angles in roll and pitch, the periods of oscillation approx- 
imate those predicted by the linearized equations of motion. The linear component of 
the yaw, dependent on the initial angular velocity, is also evident. 





Figure 3.3 Simulated Roll, Pitch, and Yaw 
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Chapter 4 


Linear Covariance Analysis 

This chapter describes the analytic development of and some of the implementa- 
tion issues involved in the linear covariance analysis of the RADCAL interferometer. 
An overview of the generic dynamics, measurement, and update equations is first pro- 
vided. Also outlined are the analytic descriptions of the error budget and suboptimal 
filtering. The derivations of the two satellite selection criteria used in this investiga- 
tion, GDOP and ADOP, are included. The error state and its dynamics as well as the 
measurement equation, specific to the gravity gradient stabilized spacecraft considered 
here, are described along with a brief discussion of their implementation. 

4.1 Overview of Filter Equations 

An underlying tenet of the linear covariance analysis is that the system dynamics 
are linear or can be linearized, and can therefore be described by the following dynam- 
ics equation: 

x = Fx + w (4.1) 

where x is the state vector, F is the dynamics matrix, and w is the process noise. In the 
analysis, we are concerned not with the full state estimate itself, but rather with the 
evolution of the second-order statistics of the errors in that estimate (i.e., its error 
covariance). These statistics can be propagated in time with the knowledge of the 
dynamics given by Equation (4.1) and updated with measurements given the linear, or 
linearized, relationship of the state to the observation, z k , at time, t k : 

= H fc :x fc + v fc (4.2) 
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where is the measurement sensitivity matrix, and v*. is the measurement noise. 
Equations (4.1) and (4.2) provide the bases from which the subsequent forms of the fil- 
ter equations are derived. 

4.1.1 Standard Filter Equations 

The true state, x, is not known, however, so the best that can be done is to maintain 
an estimate of it, x . Define the error as the difference between the true and estimated 
states: 

Ax = x-x (4.3) 

The error is assumed to be unbiased, i.e., the expected value is zero. The error covari- 
ance matrix, describing the second-order statistics of the error components and the 
correlations between them, is then: 

P = E [ (Ax) (Ax)H (4 4) 

where E[ ] is the expected value function. Given an initial covariance, P 0 , equations 
(4.1) and (4.2) can be used to compute future values of P - this is the covariance anal- 
ysis. The analysis is divided into two parts: 1) the propagation step and 2) the mea- 
surement update step. 

Due to the presence of process noise which is used to allow for any inadequacies in 
the model of the true process, the covariance will grow over time if no measurements 
are taken. The propagation of the covariance matrix P between measurements in con- 
tinuous time is described via the matrix Ricatti equation: 

P = FP + PF r +Q (4.5) 

where Q is the power spectral density matrix characterizing the process noise. This 
quantity is related to the noise covariance as follows: 

£[w(r) w(x) r ] = Q(r)5(f-T) (4.6) 
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where 8 is the Dirac delta function with units of 1/time. The noise covariance repre- 
sents the state uncertainties which accumulate over one time step due to errors in mod- 
elling the system. 

To reduce the state uncertainties, we can incorporate a measurement and update 
the state estimate as well as the error covariance (although covariance analysis is only 
concerned with the latter). We seek to use the measurement, z k , to improve the state 
estimate by a linear combination of the noisy measurement residuals applied to the 
prior estimate: 

= **+fc*[**-H*X*J (4.7) 

where the “+” indicates immediately after the update and K* is (for the moment) an 
arbitrary gain or weighting factor. The covariance of the a posteriori error. 
Ax; = x t - x k , is: 

P 4 + = E[( Ax* + ) (Ax-)n (4.8) 

Expanding this out using Equation (4.7) gives, after some manipulation, the Joseph 
form of the error covariance update: 

p; = (I - K*H*) Pjfc (I — * *H*) 7 + K*R*k[ (4.9) 

where is the measurement noise covariance. The measurement noise is assumed to 
be a white sequence with a known covariance structure. It is also assumed that there is 
no cross correlation between the measurement noise and the process noise [5, 10]: 



i = k 
i & k 


(4-10) 


E [w (?) v£] = 0 for all t and k 


(4.11) 
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We wish to find the gain factor, K*, which will give the optimal update. Defining 
this optimality as the minimum mean-square error, or the minimum trace of P, gives 
the Kalman gain, K k : 

K t = pa^h.pX+R,)" 1 < 4 '12> 

If the elements of the measurement error vector, v k , are not correlated with each other, 
the measurement noise covariance matrix then assumes a diagonal form, and each sca- 
lar measurement can be incorporated sequentially. In this case, the measurement sensi- 
tivity and weighting matrices reduce to vectors, the measurement noise covariance 
matrix reduces to a scalar, a 2 , and the Joseph form of the measurement update 
becomes: 

P k = ^1 -k^Jp^I - kjhj.J + akA (4.13) 

and the matrix inverse operation becomes a simple division: 

k, = (4.14) 

h»PA + « 

Substituting Equation (4.14) into Equation (4.13) gives the standard form of the opti- 
mal update: 

P k = ^k~ s k^k^k (415) 

where 

s k = h * p A + a 

The recursive loop for the optimal error covariance is summarized in the figure 
below. Note that the P matrix does not depend on the actual measurements themselves, 
just their statistics and the measurement geometry. Multiple measurements are pro- 
cessed sequentially using the scalar update form in Equation (4.15). 
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Po 


Multiple scalar 
updates 




Compute gain & update P: 

h = h I p A + “ 2 


k. = 


P A 


P + = P -t k k r 

r k r Jc J Jc K Jt K k 


Propagate: 

P = FP + PF r + Q 


Figure 4.1 Recursive Loop for Optimal Error Covariance 

The optimal Kalman filter gain minimizes the error covariance at the present mea- 
surement update as well as during the following extrapolation [22]. In other words, the 
gain which yields optimum performance at some later point in time is the same gain 
which yields optimum performance at the present measurement update. 

4.1.2 Error Budget 

Using consider state analysis, all of the error sources in the “truth” model, whether 
modelled or not in the onboard filter, affect the errors in the filter estimate. The error 
budget is a tabulation of the relative contributions of individual error sources, or 
groups of error sources, to system accuracy. In one form of the error budget calcula- 
tions, the time history of the filter gain is recorded and then processed repeatedly in 
determining the individual effects of each error source; multiple runs are therefore 
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required. The ensuing equations provide a more elegant method in tabulating an error 
budget with a single run [24]. 

We wish to decompose the total error covariance, P Tota i, into its constituents: 

? Total = + < 4 ' 16 ) 

The first term on the right is the contribution from the initial conditions, the middle 
term is the contribution from the error states’ process noise, and the final term is the 
contribution from the measurement uncertainty. O is the state transition matrix from 
the initial time, t Q , to the final time, tf Note that there can be as many P q and P M terms 
as desired. The quantities of interest are the effect of the above contributions on the 
final roll, pitch, and yaw mean-square errors - the first three diagonal terms in each 
component will therefore be examined. Note that in the form given in Equation (4.16), 
the initial condition contribution to the final attitude error variances is the contribution 
from the entire initial covariance, not just the initial attitude error variances. 
Equation (4.16) will require a slight modification to extract this information subset 
(see Appendix A). Continuing with the general equations of the error budget, each 
component of the covariance matrix is propagated and updated, though in a slightly 
different form from the total covariance matrix, as shown in the following sections. 

4.1.2.1 Propagation 

The usual discrete form in extrapolating the error covariance between measure- 
ments is, dropping the subscript, k, for convenience: 

P + = cpP(p r + N ( 4 . 17 ) 

where P and P + are the covariances just after the last measurement update at t k _^ and 
after extrapolation to t k , respectively, and cp = cp (t k , t k _ { ) is the state transition 
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matrix between the two times (not to be confused with the total phase measurement, (()) 
given by the relationship: 

x ('*) = 

where for a stationary process in which F is constant [5, 10]: 


The last term in Equation (4.17) is the process noise covariance over the time interval 
and is related to the spectral density, Q, by [5, 10]: 


N = f 


(f*T) dx 


(4-19) 


(See Appendix A for the details on calculating N.) Substituting Equation (4.17) into 
Equation (4.16) gives the extrapolation of each component of the covariance: 

<TP 0 <I> +7 + P e +K = (HOPo^ + Pg + P^ + N 
The extrapolation of each component is seen from the above equation to be: 


0 + = <])<D 

P q = (JjPq^+N (4.20) 


4. 1.2.2 Measurement Update 

The Joseph form of the measurement update resembles the extrapolation step 
shown in Equation (4.17): 

P + = SPS T + M ( 4 . 21 ) 

where, for the scalar update with an arbitrary weighting, k : 
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S = I-kh r and M = cx 2 kk r 

Thus, similar to propagation, the update of each component of the decomposed covari- 
ance matrix takes the following form: 

= S<D 

P* = SP q S t ( 4 . 22 ) 

K = sp^+m 

4.1.3 Subop timal Filter Equations 

The ultimate goal of covariance analysis is to design a filter for an onboard naviga- 
tor. The previously described error budget provides a breakdown of the relative contri- 
bution of each error source to the final state uncertainties and allows the states of lesser 
consequence to be identified. Due to hardware constraints, i.e., limited storage capac- 
ity and processing power, it is often necessary to reduce the order of the filter used in 
an onboard system. We therefore seek the suboptimal filter giving the best perfor- 
mance, or the minimum-variance reduced-order (MVRO) filter. To assess the perfor- 
mance of this suboptimal filter, we can cycle the suboptimal gains through the truth- 
model update equation, in which the states to be omitted are not estimated. Their 
effects, however, are felt in the states that remain; these omitted states are thus referred 
to as “consider states.” 

Equation (4.15), since it represents an optimal update, is inappropriate for rating 
the suboptimal filter’s performance. The Joseph form in Equation (4.13) is used in pro- 
cessing a suboptimal weighting vector, k, and can be reworked into a form which 
requires fewer operations and is intuitively clear [23]: 

P + = P-5kk r +5(k-k) (k-k) r ( 4 . 23 ) 
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where the subscripts have been omitted for convenience. From Equation (4.23) we 
readily see that any gain which is not optimal results in worse performance. 

In consider state analysis, the nxl optimal weighting vector, k, is computed, and 
then zeros are inserted into the components of the n-m states that will not go onboard: 

nx 1 optimal gain, k 
f K ^ 

k = WW 7 Ph(h r Ph + a 2 ) 1 
>- m elements of k 

- 
— 

_ n-m zeros for 
consider states 



where 


W = 


I 



This weighting vector is an “optimal-suboptimal” gain to be used in Equation (4.23) in 
computing the actual error statistics of the “best” reduced-order, or optimal-subopti- 
mal, filter. The term “optimal-suboptimal” is somewhat of a misnomer since this filter 
will give minimum variance in a global sense, not necessarily at every time. That is, 
another suboptimal gain can be computed to give a smaller variance at a particular 
time [22], This other gain, however, will result in performance which is worse than 
that of the optimal-suboptimal overall. The optimal-suboptimal algorithm provides the 
best performance, in general, which can be achieved in reducing the order of the filter. 
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and is useful as a benchmark for “tuning” the actual suboptimal filter to go onboard. 
The process of consider state analysis is summarized in Figure 4.2. 


Po 


Multiple scalar 
updates 



Figure 4.2 Consider State Analysis Block Diagram 


The actual onboard filter will not give the optimal-suboptimal performance 
because it does not consider the omitted states. To analyze the performance of this 
reduced-order filter, its “optimal” weighting vector, k sub , is computed and padded with 
zeros for use in the “truth” filter: 


k = Wk 


sub 
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Figure 4.3 shows the major steps in the suboptimal filter analysis. This time, however, 

two covariance matrices must be computed: 

P sub : mxm covariance matrix of the suboptimal filter. Its only relevance is in 

generating the suboptimal gain sequence. 

P nxn covariance matrix of the truth model based on suboptimal filter 

gains. It represents the statistical performance of the suboptimal filter. 




Sequence of suboptimal 
error covariances 


Figure 4.3 Steps in Suboptimal Filter Performance [5] 
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The suboptimal filter can then be “tuned” to approach the best performance given by 
the optimal-suboptimal filter. This tuning, which usually involves measurement under- 
weighting and adjusting the process noise, will most likely require multiple runs 
through the process shown in Figure 4.3. 

Performance results of optimal-suboptimal filtering using consider states are dis- 
cussed in Section 5.4. 

4.2 Satellite Selection Algorithms 

Two criteria were used to rate the effect of GPS satellite geometry on the attitude 
solution: Geometric Dilution Of Precision (GDOP) and Attitude Dilution Of Precision 
(ADOP). GDOP is primarily used for translational state estimation and is considered 
here only because the satellites that RADCAL selected were chosen based on this fig- 
ure-of-merit. Consequently, the satellites in RADCAL’s “best view” may not necessar- 
ily be best in an attitudinal sense. For this reason, the ADOP figure-of-merit was also 
employed in the linear covariance analysis for purposes of comparison. A brief expla- 
nation of the two criteria is provided herein. 

The DOP factor is a measure of how satellite geometry degrades accuracy - of user 
position and time estimates when using GDOP, or of user attitude when using ADOP. 
It is, in effect, the mapping of measurement errors into state errors due to the effect of 
satellite geometry: 

c x = DOP ■ c v (4.24) 

where v is the measurement error and x is the RSS state error. The DOP scaled by the 
measurement noise is the instantaneous RSS state error using only the measurement 
and no a priori knowledge. Satellites are chosen so as to minimize the DOP. 
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The DOP is derived as follows. Let each measurement error be zero-mean, have a 
unity la error, and be uncorrelated with errors in measurements to other satellites. 
Assuming the system is not underdetermined, Equation (4.2) can be solved for the 
error state vector, x, using the method of least squares. The error state covariance is 
then given by: 

Cov (x) = E [xx r ] = (H r Cov(v) _1 H)"‘ = (H r H)~‘ (4.25) 

since Cov(v) is an mxm identity matrix, where m is the number of satellites in the best 
set. In the above equation, H is the mxn measurement sensitivity matrix: 


H = 


^ obs , obs , .. 

d;c. 


dx, 
d £, 


3 ,, 

•=; — obs„ 
dx l 


. — obs , 

ox„ 1 


.. — obs m 
dx n m 


(4.26) 


where obsj is the i th observation or measurement and Xj is the j th component of the 
state - H will differ depending on the DOP type. For GDOP, the measurement is the 
pseudo-range, p, which can expressed as a function of user position and clock error: 

Pobs = \r GPS -Tuser\ +C &user ( 4 - 27 ) 

Row i of H in this case is merely the unit line-of-sight vector to satellite i, the sensitiv- 
ity of the pseudo-range measurement to user position, and a 1 for the clock bias term: 


H, = 


p i pi pi 1 

1 2 e 3 


(4.28) 


RADCAL’s best selection consists of a set of up to six satellites, thus making H a 6 x 4 
matrix in computing GDOP. 
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For ADOP, on the other hand, the observation becomes the single, double, or triple 
phase difference. Ignoring integer ambiguity, line bias and multipath for the moment, 
the single phase difference to satellite j is: 

A(j> ; = Ar = bV (4.29) 

Recalling Equation (1.3), the sensitivity of the phase difference to the i* component of 
the attitude is: 


- b 


r9M 


90 


BE ) 

27 


(4.30) 


and the j th row of H is: 


H ) = 


3 -Aty irAV 3 AV 


90, 90 2 


90, 


(4.31) 


With six satellites in the best view and three baselines, as in the RADCAL scenario, H 
then becomes an 18 x 3 matrix. Section 4.4.1 addresses the above measurement sensi- 
tivity at greater length. 

Having computed the measurement sensitivity matrix, we can calculate the DOP 
factor, which is defined as the square root of the trace of Cov(x) when Cov(v) is iden- 
tity: 

DOP = Jtrace [ (H r H) 
that is, for the positioning problem: 

GDOP = Jot + ct + cl + o* At 
and for the attitude problem: 

ADOP = Jol + ol + G e 2 3 (4.34) 


(4.32) 


(4.33) 
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GDOP is unitless, or [length/length]; so a pseudo-range measurement uncertainty in 
meters maps directly into a position-time uncertainty in meters. ADOP, however, maps 
the phase uncertainty into attitude uncertainty - its units are therefore slightly different. 
The phase uncertainty is given as a fraction of a cycle, or length in meters. The units of 
ADOP are then rad/m. 

If the measurements are correlated, Cov(v) is not identity, and the DOP is weighted 
accordingly. Section 4.4.2 discusses this measurement correlation.As we shall see 
later, choice of the DOP type will affect the satellites chosen and, consequently, the 
evolution of the attitude errors. 


4.3 The Error State and Its Dynamics 


The state of primary interest is, of course, the attitude error - the error in roll, pitch, 
and yaw angles from the LVLH frame. We must include, however, as many other 
states as we can think of and reasonably model which affect that attitude error. For this 
investigation, the entire error state vector is chosen to be: 


x 



60 r 5T r Sl[ ... 5lJ 5(3, ... 50* 

ant c 



( 4 . 35 ) 


where 


Ob = attitude error vector (3x 1 ) 

50 = attitude error rate vector (3x1) 

5T = torque error vector (3x1) 

51, = lever arm vector, i.e., vector from origin to i th antenna in body frame (3x1) 

5Pj = path delay error from i th antenna to receiver (lxl) 

£, = multipath vector (24x1) 
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where N ant = number of antennae = 4. The individual lever arm vectors to each 
antenna are chosen as states to avoid modelling the spatial correlations that exist 
between baseline vectors due to the common master antenna. The same reasoning 
applies to the path delays. Note that the full 1 th path delay may not be zero, but its error 
is assumed to be zero mean. Let there be one multipath state per antenna per channel. 
The Quadrex receiver on-board RADCAL has N chan = 6 channels, giving rise to 24 
multipath states: 




1 1 1 1 2 
1 2 3 4 1 




( 4 . 36 ) 


The subscript and the superscript in the above definition refer to the antenna and satel- 
lite (channel), respectively. The error state therefore consists of 49 elements total, 
whose time evolutions are described by the system dynamic equation, given in 
Equation (4.1). 

The attitude dynamics result from the linearization of Euler’s moment equations 
including the gravity gradient torque. The error torque is then due to any residual 
torques not modelled including higher-order gravity gradient, solar radiation, atmo- 
spheric drag, and other perturbative affects, characterized as a Markov process with a 
half orbital period time constant, or approximately 3000 s for RADCAL’s 815 km alti- 
tude orbit. The lever arm and path delay errors are also characterized as Markov pro- 
cesses with time constants equal to that of the error torque. Each multipath state is also 
modelled as a time-correlated Markov process, with the time constant defined as the 
time required for the LOS to the satellite to rotate through 30 degrees with respect to 
the vehicle. There are then six multipath time constants, corresponding to each of the 
six satellites being tracked. For any one satellite, however, the multipath to the four 
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antennae are spatially correlated, the degree of correlation quantified by the coeffi- 
cient, p^. 

It is convenient to partition the dynamics matrix as such: 


F = 


Fl 9 ,9 ^ 9*40 

.^40*9 ^ ^40*40] 49*49 


( 4 . 37 ) 


where FI depends on the attitude error dynamics, and F2 defines the dynamics of the 
lever arm, path delay, and multipath errors. These will be discussed separately in the 
next two sections. 


4.3.1 Attitude Error State Dynamics 

We wish to derive the form of FI, which defines the dynamics of the state sub- vec- 
tor: 


xl 


86 

5'e 

8T 


( 4 . 38 ) 


where 80 is the vector of error angular displacements about the nominal body frame 
axes, and 8T is the vector of error torques about those same axes and causes attitude 
error accelerations. Then FI is: 


FI 


3(80) 8(80) 8(88) 
8(80) 8(80) 8(8T) 
8 ( 88 ) 8 ( 88 ) 8 ( 88 ) 
8(80) 8(80) 8 (ST) 
8(8T) 8 (ST) 8(8T) 
8(80) 8(80) 8(8T) 


( 4 . 39 ) 


which results from linearizing the dynamic equations given in equations (3.3) and 
(3.4) and repeated here: 
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4«D) = I;'(T-(0X(I M C0)) 


£(M, b ) = M, s (a>x] 

The derivation of the 6x9 upper block of FI was given in Section 3.2.1. Now only the 
sensitivities of the error torque need be obtained. This will be easy since the error 
torque is modelled as a Markov process. The dynamic equation of the error torque 
then takes the form: 

5T = — 8T + w r (4.40) 

% T 

where t t is the half-orbital-period time constant and w T is the noise on the error 
torque process. The sensitivity of the error torque rate to the error torque is then sim- 
ply: 


9 (St) -4- 

3 (ST) x T 


(4.41) 


where I is the 3x3 identity matrix. 

We now have all nine 3x3 block matrices which comprise the attitude error 
dynamics matrix, FI, shown in Equation (4.39). In summary, this matrix is shown: 



®3x3 

^3x3 

®3*3 

FI = 

^■3x3 


^•3*3 


®3x3 

®3*3 

~ l 3x3 

l r 


(4.42) 


where A, B, and C are given in expressions (3.17), (3.18), and (3.19), respectively. 
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4.3.2 Lever Arm, Line Bias, and Multipath State Dynamics 

We now focus our attention on the dynamics of the second state sub-vector, which 
consists of the measurement-type states: lever arm, line bias, and the multipath errors. 
Thus, 


x2 


T 



(4.43) 


Each of these elements are modelled as a Markov process. The time rates of change of 
the lever arm and line bias errors can then be expressed as: 


d(Sl«) _ _i T g 
9(51.) ‘j 

d(5p/) _J_g 

9(5p.) t p '' 


(4.44) 


(4.45) 


where x x = tp = half-orbital-period and 5y is the Kronecker delta function which is 1 
for i = j and 0 for i * j for i,j = 1, ..., 4. Modelling the lever arm and line bias states in 
this way does not take into account the correlations resulting from the physical phe- 
nomenon assumed to drive these errors. The expansion caused by heating, for 
instance, will cause all four lever arms to expand outward in concert. So the model 
used here assumes a pessimistic performance in neglecting these correlations. 

The final set of states, multipath, are also modelled as exponentially correlated (in 
time) random variables. Since we have assigned one multipath time constant for the 
antenna array for each satellite, the time rates of change of the multipath states for the 
four antennae to satellite i are given by: 


= — — 1,_„ for i = 1, ..., 4 

3(4') V 


4x4 


(4.46) 
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The multipath time constant, is computed at each measurement time by calculating 
the angular velocity of the selected satellite relative to the body frame. 

The dynamics matrix for these states is: 


F2 = 


-*i 


12x12 


0 




4*4 


0 


-±I 


4*4 


0 


0 


-±i 

V 


4*4 


( 4 . 47 ) 


40*40 


We now have the constituents of the entire (49x49) dynamics matrix, F. 


4.3.3 Process Noise 

The process noise vector, w, accounts for inadequate modelling of the system 
dynamics. To be on the conservative side, a zero-mean Gaussian white noise process is 
incorporated into the attitude acceleration process, resulting in a random walk compo- 
nent in the attitude rate. The nominal value of the noise power spectral density for the 
three components of the acceleration is 10' 7 (deg/s 2 ) 2 /Hz, resulting in a random walk 
in attitude rate with a standard deviation growth of 3.2xl0‘ 4 (deg/s)/sqrt(s). This level 
of noise is chosen to give steady-state rms attitude error rates in the nominal run that 
are commensurate with the results given in Reference [2], 

The torque, lever arm, and line bias error states are modelled as Markov processes, 
or exponentially correlated random variables (ECRV’s). The white noise spectral den- 
sity, which generates the process described by the autocorrelation function of the 
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ECRV in Equation (2.1), is given in terms in terms of the variance of the state, x, as 

[ 10 ]: 



The multipath error states are also assumed to behave as Markov processes, though 
their noise terms differ from those of the torque, lever arm, and line bias error states in 
that spatial correlations between the antennae are assumed to exist. To determine the 
form that the multipath white noise PSD will take, let us consider the multipath to two 
antennae from one satellite and prescribe that the steady-state covariance of these two 
states is of the form: 


where 


2 2 
2 2 


= multipath correlation coefficient 

Insert this into the matrix Riccati equation with unknown dynamics and process noise 
matrices: 


P = FP + PF r +Q 


00 = 2 f a f b °\ P^\ + q a q b 
9 0 . fb fa p^G^ G^ Qb Qa_ 

Now if it assumed that each multipath state has the dynamics of a typical ECRV, both 
having the same time constant, then/ fl = -l/x^ and f b = 0, and the above equation can 
be used as a constraint to solve for q a and q b : 
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The observation or measurement of the process occurs at discrete points in time as 
dictated by the linear relationship: 

z * = H * x * + V * 

where, for the case considered here, 
z k = measurement vector (3x1) 

H k = measurement sensitivity matrix (3x49) 
v k = measurement noise vector (3x1) 


The measurement vector used is the difference between the observed and predicted 
phase differences: 
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(4.49) 



A<t>! 


A4>, 

II 

N 

a<> 2 

- 

A<i>2 


A<h 


_A<j>3 


The predicted values are denoted with a “ A ”, and the subscript refers to the baseline. 
There will be up to six of these measurement vectors per epoch, one for each satellite 
in view. Note that it should be possible to take a measurement if any two antennae 
view a particular satellite, but for some reason RADCAL only processes measure- 
ments to a satellite if it is in view of all four antennae at once. Nonetheless, the system 
is overdetermined with 18 phase difference measurements per epoch. 

The actual phase difference between antennae 1 (master) and i (slave) to satellite j 
can be expressed as: 

=A +p 1 -p. + ^ I -^ + v , 1 -v' (4.50) 

Equation (4.50) will be used in the following section to derive the sensitivity of the 
phase difference to each of the error states, with the integer ambiguity, n, assumed to 
be known. The final two terms in the expression, which comprise the measurement 
noise, will be addressed as well. 


4.4.1 Measurement Sensitivity 

Let us consider just one phase difference so we need only deal with a measurement 
sensitivity vector: 


T 3 (A4>) 
3x 


a(A<J>) 3(A<|>) 3(A<(>) 9 ( A<t>) 3(A<t>) 3(A<{>) 
96 30 3T 31 9(3 31; 


or 


h 


T 



h- h r h ' h l h 


r 

§ 
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where it is understood that there are four terms each for the lever arm and line bias cor- 


responding to each of the four antennae. We proceed to derive each term in the sensi- 
tivity vector. 

Recalling the range difference from Equation (1.3): 

A r = b B M B£ e £ (4.51) 

where b is the baseline vector in body coordinates, e is the LOS in ECEF coordinates, 
and M 5 £ is the transformation between the two frames in which is embedded the roll, 
pitch, and yaw from the LVLH frame. Express this true transformation in terms of a 
transformation to the nominal frame followed by a correction transformation: 

M ,e = M bb M 6s = <>- [SexDM^ 

and substitute it into Equation (4.51) to get: 

Ar + 5 (Ar) = (b B + 5b B ) T M hE e E -bl [80x] M~ £ e £ 

The perturbation in the range difference is therefore: 

8 (Ar) = 5bjM^ £ e £ + b B [ (M^ £ e £ )x]50 

From this expression, we can get the sensitivity to both the attitude and the lever arms 
in use: 


,r _ 3(A<|>) _ 3 (A r) _ r 

“e ~ -m _ ^ — t 


39 - 30 - ^(Mfl £ e £ )x]b B = ~[(M- £ e £ )x] (1,-1,) 


h ' r = ^ = tt 1 = -*x* and k = “-ir = ^ 


t _ a (A<|>) _ a (Ar) _ 

- E b * 

‘1 UM 1 

where we have used the fact that the baseline of interest can be expressed in terms of 
the nominal and error lever arms in body coordinates: 


b B = a + 5i,)-(i, + 8i.) 

Referring once again to Equation (4.50), we see that for the line biases: 
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h _ a(A4>) _ J d h _ 3(A<t>) _ _J 

" ap, " p- " ap, " 1 


and for the multipath: 


h . = <^4> = i and h, = = -1 

« 34 ', « 34 ' 


and the sensitivities associated with the other antennae and other satellites are zero. 
The full 24 element multipath term for antennae 1 (master) and 2 (slave) to satellite 2, 
for example, is: 


^ - [®ijc 4 1-10 0 0 1_(4 0 Ix4 0 lx4 Oj 

Finally, we note that the phase difference is not sensitive to the attitude rate nor the 
torque: 


t _ a(A(j>) _ 0 d h r _ a(A(|)) _ 0 
n. - ^ - u u3 ana n T - ^ - u 


1x3 


We now have our fully populated 49x1 measurement sensitivity vector. 


4.4.2 Measurement Noise 

The measurement noise is assumed to be a white sequence with a known covari- 
ance structure. It is also assumed that there is no cross correlation between the mea- 
surement noise and the process noise [5, 10]: 


ECv.vD = ( R ‘ ‘ = k 

L * ‘ VO f * k 


(4.52) 


£[w(0v 3 = 0 for all t and k (4.53) 

The noise term in Equation (4.50) for the phase difference is a combination of the 
total phase measurement noises from the master and slave antennae. The noises on the 
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three phase differences will be correlated due to the common master antenna. Let us 
first make some assumptions about the nature of the total phase measurement noise: 

£[ v N,*] = ( = 0 for i *j ) 

£[vi v *] = cl 6* ( = 0 for j*k) 

<> j 

The first indicates that the measurement noises of two antennae to the same satellite 
are independent, and the second states that the noises between channels are indepen- 
dent. Based on the first assumption, we can say that the 4x4 total phase measurement 
noise covariance associated with each satellite, R^, is diagonal: 

R * = ^[VD = * 1 } 4,4 

To characterize the correlated measurement noises of the phase differences, we can 
use the single difference operator, SD, defined as [12]: 


A<j)~ 


-110 0 

Atj^ 

= 

-10 10 

A< t>3 


-1 0 0 1 


4>i 

<t>2 

(J ) 3 

4>4 


SD<]> 


Then the phase difference measurement noise vector becomes: 


V A$ = AV $ = SDv * 

and the 3x3 measurement noise covariance of the differences is: 

= E[v M <J = ftSDv^JSD 1 ] = SDR,SD r 

This operation results in a covariance matrix with off-diagonal terms: 


R ai = o;SDSD r = 



2 1 1 
1 2 1 
1 1 2 
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nr 

(The matrix, SD SD , is used in weighting the ADOP, as discussed in Section 4.2.) The 
sequential filter that we wish to use to process one phase difference at a time, however, 
requires uncorrelated measurements. We therefore need to decorrelate the measure- 
ments by diagonalizing R A()) , as suggested by Sullivan [26]. 

Any nxn symmetric matrix. A, can be diagonalized into a matrix, D, by construct- 
ing a matrix, U, the columns of which are the eigenvectors of A and are linearly inde- 
pendent. The non-zero elements of D are the eigenvalues of A. The matrix A can then 
be diagonalized as follows [25]: 


U~*AU = D 



0 


0 


where the ^’s are the eigenvalues. Unit eigenvectors give an orthonormal matrix U; 
therefore, U" 1 = U r . Using this logic, we can diagonalize R A()> : 


U r R A ,U = 


a] 0 0 

0 04 0 

0 0 aj 
_ 


Decorrelating the phase difference measurement noises in this way has altered the 
measurement equation. Since 


U r R A „lJ = E[u r v A *vI*u], 

the new vector of uncorrelated phase difference measurements is: 

z' = U 7 z = U r Hx + U 7 v A * 

and each scalar measurement can now be processed sequentially by realizing that the 

T T 

new measurement sensitivity vectors, h' , are the rows of U H. The 3x3 modal 
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matrix, U, can be precomputed since the measurement statistics do not change. Pro- 
cessing scalar measurements, as opposed to incorporating the three correlated phase 
differences in batch, results in about a 50% reduction in computation time for six-sec- 
ond propagation and measurement intervals. 

4.5 Simulation Overview and Implementation Issues 

The flow diagram shown in Figure 4.4 illustrates the major steps in the covariance 
analysis. First, the initial conditions are set. These include definition of the nominal 
trajectory, which is a circular orbit and a local vertical attitude hold. The spacecraft is 
nominally in a polar orbit and starts at zero degrees latitude. Also initialized are the 
error covariance, process noise spectral density, and measurement noise covariance 
matrices. 

The GPS constellation and the spacecraft states are then propagated to the current 
time. The GPS satellite positions are found by interpolating over the precise ephemer- 
ides made available by the International GPS Service, rather than from a simulated 
GPS constellation. The precise ephemerides files contain the ITRF-93 positions, accu- 
rate to within 10 cm, of all active satellites at regular time intervals. This high degree 
of position accuracy is not required for attitude determination (e.g., only a three degree 
error in phase difference, about 1.6 mm, over a one meter baseline results from a posi- 
tion error of 30 km). The precise ephemerides are used in anticipation of filtering the 
actual flight data and are easier to manipulate than the broadcast ephemerides in 
RINEX format. The file obtained for this investigation corresponds to GPS week 708, 
day 6, or August 7, 1993, that of the RADCAL flight data. (See Appendix B for dis- 
cussion of RADCAL data.) Calculating the user satellite state is trivial for the circular 
orbit and LVLH attitude. 
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If it is time for a measurement, the best group of satellites are selected based on 
either GDOP or ADOP, described earlier. GDOP is the default criterion since RAD- 
CAL’s best view is defined in this way. The multipath states in both the covariance and 
process noise spectral density matrices are then updated. If a satellite in view has 
switched or dropped out since the last measurement time, the multipath elements cor- 
responding to that channel must be reinitialized. The 4x4 block diagonal matrix in the 
covariance matrix is set to its initial, steady-state value as given in Section 4.3.3, and 
the correlations with all other states are zeroed. If the satellite was replaced, the new 
multipath time constant is calculated from the relative angular velocity. If the satellite 
dropped out and was not replaced, the multipath variances for that channel will remain 
at their steady-state values, since no measurements are taken. The noise spectral den- 
sity associated with that channel must also be updated, since it depends on the time 
constant. 

With the satellites selected and the multipath states updated, the measurement sen- 
sitivity and thus the weighting, optimal or otherwise, can be calculated. The measure- 
ment can then be incorporated into the total covariance matrix via Equation (4.23). 

The dynamics matrix given in Section 4.3 is computed, and the covariance is then 
propagated forward to the next time. The propagation of the total error covariance 
matrix as described in Section 4.1.1 is accomplished in this investigation via fourth- 
order Runge-Kutta integration of the Ricatti equation. The dynamics matrix is 
assumed to be constant over the integration time - a small time step compared to the 
system dynamics is thus required. Upon propagation, the process repeats until the des- 
ignated final time. 

Some of the implementation issues involved in performing the error budget are 
discussed in Appendix A. 
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Chapter 5 


Simulation Results 

The following sections present the results obtained from the simulation outlined in 
the previous chapter. A nominal run is defined and its results are discussed. The per- 
formance sensitivity to various parameters such as orbital inclination, satellite selec- 
tion method, and noise levels are examined, as well. Finally, the performances of two 
optimal-suboptimal filters are discussed. 

5.1 Nominal Run Description 

Table 5.1 shows the nominal run assumptions consistent with the RADCAL sys- 
tem. The error terms are, of course, only a best estimate and will be varied in a later 
section to determine performance sensitivity to these assumptions. States with associ- 
ated time constants (x x ) are modelled as Markov processes. The performance results 
presented in this section are fully optimal and will show what might be possible 
onboard given enough computation power. They will provide insight into system 
behavior. 


Table 5.1 Nominal Run Conditions 


Parameter 

Condition 

Orbital Altitude (circular) 

815 km 

Orbital Inclination 

90.0° 

Attitude Hold 

LVLH 

Antenna Baseline 
Configuration 

3 Baselines as described in 
Section 1.4 

Antenna Half-Cone Angle 

85.0° 

Antenna Boresite Cant Angle 

17.5° 

Satellite Selection Criterion 

GDOP 
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Table 5.1 Nominal Run Conditions 


Parameter 

Condition 

Principal Moments of 
Inertia 

26.400, 26.400, 5.813 kg m 2 

Initial Attitude, Attitude Rate 
Uncertainty (each axis) 

5.0°, 0.1°/s 

Torque Error 

ax = 5.0 xl0‘ 7 Nm, 
Tx = 3035 s 

Lever Arm Error 

aj = 2 mm, Tj = 3035 s 

Line Bias Error 

ap = 2 mm, tp = 3035 s 

Multipath 

a^ = 5 mm, » 500 s 

Multipath Correlation 

50% 

Attitude Error Rate Process 
Noise PSD 

10' 7 (°/s) 2 /s 

Total Phase Error 

CO 

o 

o 

Propagation Step Size & Time 
Between Measurements 

6 s 


The nominal conditions in the orbital altitude and inclination, measurement fre- 
quency, antenna baseline configuration, antenna boresite cant angle, satellite selection 
criterion, and principal moments of inertia are all chosen to be consistent with those of 
RADCAL. The 85° antenna half-cone angle chosen for the baseline case is based on 
GPS satellite visibility. A satellite is considered in view if it is visible to all four anten- 
nae simultaneously. Figure 5.1 reveals the simulated visibility for the 85° half-cone 
angle to be comparable to its real counterpart. See Figure B.2 in appendix for the num- 
ber of satellites in RADCAL’s best view. At least four satellites will be in view more 
than 90% of the time in both cases.The simulated spacecraft is over the poles at t = 
1500 s and 4500 s, respectively. 

The nominal attitude is chosen to be the ideal LVLH orientation of a gravity gradi- 
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Figure 5.1 Simulated GPS Satellite Visibility: 85° Half-Cone Angle 


that the craft can be off by 20 degrees in roll and pitch, and that the yaw grows secu- 
larly in time [2, 8, 28]. The LVLH assumption will do for the purposes of linear cova- 
riance analysis. (Not to imply that the performance is insensitive to attitude.) The 
assumed nominal torque error corresponds to 5% of the modelled gravity gradient 
torque arising from 5° roll and pitch angles using the inertias shown. 

The relative antenna movement and changes in electrical path length, which intro- 
duce errors in the attitude solution, are attributed to temperature variations. In assign- 
ing values for the lever arm and line bias error processes, the following simplistic 
argument is made. The change in length of the lever arm, A /, is a function of the tem- 
perature change, AT, and the material’s thermal expansion coefficient, a: 

Al = la (AT) 

where l is the nominal length. The temperature of a typical spacecraft ranges from 
-34°C to 71°C [31]. Assuming a = 20E-6 m/(m 0 C), representative of aluminum, this 
temperature difference results in a change in length of two parts per thousand, or 2 mm 
for a 1 m lever arm. This conservative number is used for the lever arms here, which 
are less than a meter. The steady-state line bias rms error is also assumed to be 2 mm 
in the nominal case. 
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The characteristic multipath level is assumed to be 5 mm, as suggested by Cohen, 
Lightsey, and Parkinson [8, 18], and its correlation between antennae is nominally 
assumed to be 50%. 

Though less than 1 mm of phase measurement error has been suggested for RAD- 
CAL’s patch antennae [8], a slightly more conservative figure of 3° or about 1.6 mm is 
used, instead, in an attempt to absorb the error induced by the antenna’s asymmetric 
reception pattern. 

The step size shown is the average time between RADCAL measurements. In gen- 
eral, an appropriate integration step size is 1/100* of the minimum period of the sys- 
tem dynamics or smaller. Given the slow dynamics of the gravity gradient stabilized 
spacecraft as seen in Chapter 3, it would seem reasonable to assume that the attitude 
error dynamics can be considered constant over an integration time step up to 30 - 40 
s. The torque, lever arm, and line bias error dynamics behave on a similar time scale. 
The multipath, however, as modelled here, operates on a smaller time scale (higher 
dynamics) with a time constant on the order of 500 s. Results from runs with various 
step sizes, not shown here, reveal no significant difference (0.01% maximum differ- 
ence in variances over 10000 s) in performance between time steps from 0.1 to 6 s. 
The dynamics matrix is therefore safely assumed to be constant over the nominal 6 s 
time step. This is also the maximum time step possible because measurements nor- 
mally occur at 6 s intervals. 

5.2 Nominal Run Results 

5.2.1 Attitude Error and Error Rate 

Figure 5.2 shows the post-measurement rms errors in roll, pitch, and yaw and their 
root-sum-squares (RSS) for a 10000 s interval, more than an orbital period. The stan- 
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Time (s) 



dard deviations in the roll, pitch, and yaw jump quickly from their initial 5 degrees to 
nearly their “steady-state” values, limited by process and measurement noises. 

Notice from Figure 5.2 that the yaw angle is the most observable (i.e., its error is 
smallest), whereas the roll and pitch angles are less observable and have very similar 
error characteristics. This can be explained geometrically as a result of the fact that the 
three baselines lie in the horizontal plane. With this geometry, satellites directly over- 
head are best for determining roll and pitch, while satellites in the horizontal plane are 
best for determining yaw (disregarding antenna limitations and susceptibility to multi- 
path, of course). 
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Clearly, a conflict exists in choosing the best satellites for yaw determination ver- 
sus those for roll and pitch determination. By means of ADOP it is possible to resolve 
this conflict in an optimal way. Nevertheless, roll and pitch estimates tend to suffer rel- 
ative to yaw, as illustrated in Figure 5.2, simply because there are less overhead satel- 
lite opportunities. 

Figure 5.3 shows a sky map for the 10000 s of the nominal run. The arcs are the 
paths of the GPS satellites in view when looking to the zenith. Vehicle motion is in the 
+X direction, which is at 0 degrees azimuth, so the satellites move from right to left 
across the map. The +y direction is at 90 deg, and the +z direction (nadir) comes out of 
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Figure 5.3 Satellite Motion (Azimuth & Angle from Zenith) Sky Map 

the page. The radial component on the plot is the angle from zenith. Note that for this 
time interval in a polar orbit, no satellites can be seen directly overhead. The visibility 
“cone” is seen from the plot as being somewhat squared-off as a result of the out- 
wardly canted antennae. 
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The ragged appearance of the plots in Figure 5.2 indicate that the attitude errors 
may be at the mercy of the satellite geometry. To corroborate this possibility, two mea- 
sures of satellite geometry, GDOP and ADOP, are shown in Figure 5.4. Two curves are 




Time (s) 

Figure 5.4 GDOP and ADOP for Nominal Run 

shown in the ADOP chart - one is the ADOP calculated for the geometry of the satel- 
lites selected using the GDOP criterion, and the other is the ADOP using the ADOP 
criterion to select the best satellites. Note the units of ADOP: to get the instantaneous 
attitude uncertainty, the phase uncertainty is first converted to meters then scaled by 
the ADOP to give radians, then converted to degrees. 

A glance at Figure 5.1 reveals that between 800 and 1500 s, more than six satellites 
are in view. Looking to the insert in the above figure at the same time interval, we see 
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that use of the two criteria result, not surprisingly, in different combinations of the 
“best” six satellites, as evidenced by the discrepancy in the ADOP’s. The implications 
of the two figures-of-merit will be discussed further in Section 5.3.2. A comparison of 
Figure 5.2 and Figure 5.4 shows that the attitude errors tend toward the ADOP, which 
represents the “steady-state” condition at any given time, dependent upon instanta- 
neous geometry. 

Since the evolution of the attitude errors depend largely on that of the satellite 
geometry, perhaps it is more appropriate to cite the averages, minima, and maxima of 
the “steady-state” rms errors over, say, an orbital period, rather than the final uncer- 
tainties. Table 5.2 shows these values over the final 6000 s of the nominal run. The 
maximum values in each component do not occur at the same epoch. For each compo- 
nent, the average error as given here is actually the square root of the average of the 
post-measurement variances over the 6000 s. 


Table 5.2 Attitude Errors Over Final 6000 s 


Component 

Min RMS 
Error 
(deg) 

Max RMS 
Error 
(deg) 

Average 
RMS Error 
(deg) 

Roll 

0.457 

0.536 

0.496 

Pitch 

0.457 

0.533 

0.495 

Yaw 

0.271 

0.415 

0.317 


The attitude error rates and the magnitude, shown in Figure 5.5, exhibit similar 
behaviors to the attitude errors in that the yaw component is generally smaller than the 
roll and pitch components, which are very much alike. Note that the attitude rate is not 
directly observable from the phase differences, but rather is inferred from the attitude 
itself. 
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Figure 5.5 Attitude Error Rates of the Nominal Run 


Closer Inspection of Attitude RMS Errors 

Let us examine the evolution of the errors in greater detail. In Figure 5.6 are close- 
ups of the attitude rms error, its magnitude, and the ADOP for the selected satellites 
for the first 1000s of the nominal case. The Roll, Pitch, and Yaw Dilutions Of Preci- 
sion, are also shown in the figure, indicated as RDOP, PDOP (not to be confused with 
Position Dilution Of Precision), and YDOP, respectively. These component-wise 
DOP’s are simply the square roots of the diagonal components of (I! 1 !!)' 1 , i.e., their 
RSS is the ADOP. They represent the rms errors in roll, pitch, and yaw that would be 
obtained if only the measurements were used (no a priori knowledge) and the mea- 
surement uncertainty was unity. Inspection of the six events labeled on the plot will 
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provide more insight into the behavior of the attitude errors. Table 5.3 is a list of the 
descriptions of the satellite changes causing the events. 


Table 5.3 Description of Events Shown in Figure 5.6 


Event 

Number 

Event Description 

1 

Number of satellites in view increases from 4 to 5 

2 

Number of satellites in view decreases from 5 to 4 

3 

Satellite switches - 4 satellites still in view, better GDOP, worse ADOP 

4 

Number of satellites in view increases from 4 to 5 

5 

Number of satellites in view increases from 5 to 6 

6 

Satellite switches - 6 satellites still in view, better GDOP, worse ADOP 


Obviously, the addition of another usable satellite (usable in that we are not 
required to throw another satellite away) brings only a better knowledge of each state, 
as seen in events 1, 4, and 5. One would expect attitude accuracy to suffer when a sat- 
ellite drops out of view (event 2), as indicated by the ADOP. The yaw error does in 
fact start to increase, albeit slowly because of filter memory, as it tries to track the 
YDOP, but the roll and pitch errors continue to decrease. The satellite which dropped 
out of view was on the horizon, and as previously mentioned, the phase differences to 
this lower-elevation satellite are fairly insensitive to roll and pitch, but very sensitive 
to yaw; the roll and pitch DOP’s are not affected as strongly by the loss and continue 
to improve after it. The behavior of the roll and pitch rms errors reflects this insensitiv- 
ity. 

Some unexpected results arise from a satellite switch as in events 3 and 6. In each 
case, the GDOP of the new combination is smaller than that of the old combination, 
but the ADOP is larger. Each component of the ADOP is larger in event 3 as well, and 
yet the attitude uncertainties are seen to decrease sharply, more noticeably in yaw. 
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Figure 5.6 Attitude Errors and ADOP for 1000s of Nominal Run 
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These decreases at the switches are due not to the filter taking advantage of the new 
geometry with the memory of the previous geometry, but rather the correlation, or 
more appropriately, the decorrelation of the attitude with the multipath state at the 
switch. Up to the time of the switch, correlations build between the attitude error and 
the multipath - the uncertainties in each state will decrease to a point, but the two 
states can not be separated because of these correlations. When a satellite switch 
occurs, a reinitialization takes place - the 4x4 matrix in the covariance corresponding 
to the multipath from the new satellite to the four antennae is set to its initial, steady 
state, and the correlations between those multipath states and all other states are 
zeroed. The effect of this decorrelation is to temporarily add new information to the 
system, and the immediate result is a smaller uncertainty in both states. Correlations 
then build up again, and the uncertainties again tend toward a level consistent with the 
ADOP. The reason we do not see a subsequent increase is that the noise driving the 
attitude process is too small to be seen over short time intervals. One might contend 
with the realism of the model and could argue that a new satellite should not come in 
uncorrelated with the system, but this was the assumption made. 

The nature of the interaction between the attitude and multipath states can be sub- 
stantiated by examining some more cases. A larger amount of process noise feeding 
into the attitude causes the filter to weight the measurements more heavily, and as the 
noise level increases, the error magnitude approaches the ADOP scaled by the mea- 
surement uncertainty. Figure 5.7 shows the attitude errors resulting from a larger atti- 
tude rate process noise spectral density, 10' 3 (deg/s) 2 /s. The dip associated with the 
switch at event 3 is still visible in yaw, but the errors are tending more toward the 
ADOP than in the nominal run case with a smaller rate noise. A satellite switch also 
occurs at event 6, but note that the errors only slightly increase in roll and pitch. In 
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yaw, however, the uncertainty remains at the lower level, consistent with the lower 
YDOP after the switch. The reason that the decorrelation apparently has a greater 
impact on the yaw component has yet to be determined. 

Figure 5.8 shows the effect of eliminating multipath from the state by zeroing its 
uncertainty. We no longer see the sudden decrease at the satellite switch. The yaw 
error is seen to increase with the switch at 3, while the roll and pitch errors decrease. 
Observing the trend in the DOP’s in the previous discussion and since the multipath 
effect has been removed, these phenomena must be due to the lever arm and line bias 
states. 
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We conclude that the simulation results concerning the attitude-multipath correla- 
tion are consistent with the assumptions made in the model, and that the unrealistic 
phenomena shown in the plots are not due to a programming flaw. The “problem” lies 
with the model itself, which assumes that: 1) multipath is a zero-mean process, and 
2) the multipath associated with a newly introduced satellite is not correlated with the 
attitude. Taken over an extended period of time, perhaps it is zero-mean in reality, but 
at any one point, multipath is a bias. This bias is entirely dependent on the local reflec- 
tive environment and the signal’s direction of incidence with respect to that environ- 
ment. In other words, real-world multipath associated with a newly introduced satellite 
is in fact correlated with the attitude, contrary to the model’s assumption. 


5.2.2 Torque Errors 


Recall from Section 4.4.1 that the phase difference measurement is insensitive to 
the torque, as well as to the attitude rate. Figure 5.9 below illustrates this insensitivity. 
The GPS measurements do not measure torque directly - it must be inferred from cor- 
relations that have built up. Since the torques are small, the correlations are small, and 
the observability is low. After 10000 s the torque rms error magnitude decreases by 
less than 5%. But again, like the attitude error and error rate, the yaw component is 
slightly more observable than the torques about the roll and pitch axes. 


x io ' 7 




Figure 5.9 Torque Errors of the Nominal Run 


5.2.3 Lever Arm and Baseline Errors 

The filter has a difficult time distinguishing one lever arm from another given the 
phase difference measurement between the antennae. Figure 5.10 illustrates this diffi- 
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culty by showing the component errors for each of the four lever arms, denoted LAI, 
LA2, LA3, and LA4. The z components of all lever arm errors are equally difficult to 
observe - note that they are all the same and larger than the x and y components. This 
outcome is consistent with the results of the attitude, attitude rate, and torque error 
states. All four antennae nominally lie in the body x-y plane (shown in Figure 1.2), 
which nominally coincides with the local-horizontal plane. A flexure in the z direction 
for lever arms 1 and 3, which are aligned respectively with the -x and +x directions, is 





Figure 5.10 Nominal Run Errors of the Four Lever Arms 


equivalent to a rotation about the y axis, or a pitch. Similarly, a flexure in the z direc- 
tion for lever arms 2 and 4, which are aligned with the -y and -i-y directions, is tanta- 
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mount to a rotation about the x axis, or a roll. These apparent rolls and pitches are not 
easily estimated for the reason given in Section 5.2.1 and are difficult to distinguish 
from actual vehicle roll and pitch. 

Another item of note is that the component along the lever arm axis in each case 
(e.g., the x component of LAI) has a smaller rms error than the other two components. 
The phase difference is insensitive to changes in baseline length when the LOS is per- 
pendicular to that baseline. There is a one-to-one correspondence between a change in 
length and a change in phase difference when the LOS is coincident with the baseline 
axis. The plots indicate that the lever arm length is more observable than its misalign- 
ment, which is reasonable since misalignment of the lever arm can be misinterpreted 
as vehicle rotation. 

We are more interested, however, in the errors over the baselines, not the lever 
arms. Estimating the lever arm to each antenna was merely a convenience to avoid 
modelling the correlations between baseline flexures due to the common reference. 
The baseline 1 error is the error in lever arm 2 with respect to the error in lever arm 1 : 

b, = l 2 -lj=>5bj = 5I 2 — 51, 

The baseline error covariance is, by definition, the expected value of the outer product 
of the baseline errors: 

P b] = EtSb.SbJ] = £[(81,-81,) (81,-Sl,) 7 ] 

Expanding the right side in terms of the lever arm error covariances gives: 

P b] = £[61j5lJ] +£[51 2 8l3 -E[81,5l£] 

The above formulation gives the error covariance in the body frame, but it may be 
more illustrative to transform to a local baseline frame - one in which the coordinate 
axes are along and perpendicular to the baseline axis. Since the three baselines are in 
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the x-y plane, the transformation is simply about the z axis through the angle, a, from 
the body frame (subscript “ B ”) to the local baseline frame (subscript “A”): 

^ = Tib 

where 

T.b = T z (a) 

The covariance in the A frame is then: 

P,, = E[8b A 8b3 = ECT^Sb.Sb^] = T, s P b T SA 
For baselines 1, 2, and 3, a = -45, 0, +45 degrees, respectively. The square roots of the 
diagonal elements of the resulting three baseline error covariance matrices are shown 
in Figure 5.11. 

Note that for each of the three baselines, the uncertainty is smallest in the “along- 
track” component, that is, the component along the axis of the respective nominal 
baseline. The “cross-track,” or transverse, component in each baseline, in the direction 
perpendicular to the baseline in the horizontal plane, is somewhat less observable. 
These y A components, as shown on the plots, are difficult to distinguish from yaw 
motion. Notice the one-sided saw-tooth characters of the along- and cross-track com- 
ponents, similar to the profile of the yaw error itself - a reasonable result since these 
two components are in the plane defined by yaw rotation. The z A component, which is 
in the nadir direction, is seen to be the least observable for all three baselines. This 
component in baselines 1 and 3 can be construed as a combination of roll and pitch 
motions, while for baseline 2, this vertical component appears indistinguishable from 
vehicle pitch. 
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x A Along-Track Axis 

y A Cross-Track Axis 

z A = z B Z Body Axis 

k Z A 


► x B 



Note: x A and y A are different for each 
baseline. z A is nadir for all three. 


Figure 5.11 Nominal Run Baseline Errors 


5.2.4 Line Bias Errors 

The line bias errors associated with each antennae are indistinguishable, similar to 
the lever arm errors. Figure 5.12 shows that the line bias rms error for antenna 1 hangs 
up at a la value of approximately 1.6 mm. This profile is identical to that of antenna 3, 
since those two antennae lie along the same body axis. For the same reason, the errors 
for antennae 2 and 4 are identical and are both very close to those of antenna 1 shown 
in the figure. 
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Time (s) 

Figure 5.12 Nominal Run Line Bias Error for Antenna 1 

The line bias errors over the baselines can be obtained in the same manner as the 
baseline errors in the previous section - by differencing the appropriate antenna line 
biases and accounting for their correlations. These differential path delay errors are 
shown in Figure 5. 13 for baselines 1 and 2 - the baseline 3 results are the same as those 
of baseline 1 due to the symmetry of the two baselines about the x axis. The line bias 
errors for the first baseline appear to reach a steady-state value of approximately 2 
mm. The errors for the second baseline are seen to level off near 2.2 mm. The reason 
for this slightly higher value is not known. 




Figure 5.13 Nominal Run Line Bias Errors for Baselines 1 & 2 
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5.2.5 Multipath 


The multipath errors from one satellite are modelled as correlated between anten- 
nae, but the errors between satellites is assumed independent. Figure 5.14 shows the 
multipath for antennae 1 and 3 on channel 2 to be estimable to no better than about 4.2 
mm, even with the best geometry. The multipath errors on these two antennae, which 



Figure 5.14 Nominal Run Multipath Errors for Antennae 1 & 3 on Channel 2 



Figure 5.15 Nominal Run Multipath Errors for Antennae 2 & 4 on Channel 2 

lie on the x axis, are indistinguishable, as are the errors on antennae 2 and 4 on the y 
axis. The multipath uncertainties associated with the latter antennae for the same chan- 
nel are shown in Figure 5.15. The jumps to the 5 mm uncertainty level at 1600, 2300, 
5300, and 8000 s occur when there are less than 6 satellites in view. At those times, no 
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satellite is visible on channel 2, and the la multipath error defaults to its initial value. 
The nature of the first-order Markov process keeps the uncertainty at this steady-state 
value, since no measurements are taken to update it. The sudden dip near 300 s corre- 
sponds to event 3 in Section 5.2.1, in which a satellite switch occurred. As previously 
described, the decrease in uncertainty is due to the fact that the new satellite multipath 
errors are initially uncorrelated with the rest of the system. Note that the error did not 
jump to 5 mm because the errors shown are after the measurement. 

5.2.6 Error Budget 

An error budget was tabulated for the nominal run. Table 5.4 provides the percent- 
age breakdown of each component to the roll, pitch, and yaw error variances, as well 
as the absolute attitude errors in degrees after 10000 seconds. Note that the first contri- 
bution listed is the effect from the initial conditions (IC) of all states, not just from 
those of the attitude uncertainties. The fact that the percentages do not total to 100 is 
due partly to round-off error, and partly to the truncation of the Taylor series in the 
process noise covariance, N. The values quoted are after the measurement since we are 
interested in the best filter performance. If we had merely propagated to the final time 
and not taken a final measurement, the process and measurement noise terms would 
have commanded a greater percentage of the attitude uncertainties. 

Immediately apparent from the tabulation is that the contributions to the yaw are 
quite distinct from their counterparts in roll and pitch, consistent with previously pre- 
sented results. It should not be surprising that the effects of initial conditions diminish 
significantly over time. From the table, we see that the impact of the initial conditions 
on roll and pitch uncertainties after 10000 s is a fairly insignificant 2%. The yaw com- 
ponent, on the other hand, seems to be considerably more sensitive to the IC’s, which 
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Table 5.4 Error Budget for Nominal Run 



Per Cent of Total 
Variance 

RMS Errors 
(deg) 

Contribution 

Roll 

Pitch 

Yaw 

Ron 

Pitch 

Yaw 

1) Initial Conditions 

1.97 

1.94 

15.09 

0.065 

0.064 

0.115 

2) Process Noise: 



Attitude Rate 

1.34 

0.59 

1.55 

0.053 

0.036 

0.037 

Torque 

0.03 

0.00 

0.08 

0.008 

0.001 

0.008 

Lever Arm 

31.22 

31.60 

39.60 

0.257 

0.259 

0.187 

Line Bias 

23.16 

23.22 

1.01 

0.221 

0.222 

0.030 

Multipath 

39.77 

39.99 

38.49 

0.290 

0.292 

0.184 

3) Measurement 
Noise 

2.11 

2.29 

4.11 

0.067 

0.070 

0.060 

Totals 

99.60 

99.63 

99.93 

0.460 

0.461 

0.297 


constitute 15% of the total mean square error. This is partly due to the fact that the 
overall yaw error is smaller; therefore, the percentage effect of the initial conditions is 
larger. But this is not the whole story since the absolute contribution of the IC’s to yaw 
is larger (see last column, first row). This is a suspicious result since the yaw angle is 
more observable (as indicated by the smaller total uncertainty) than the roll and pitch 
angles and should therefore be less influenced by the initial conditions. To investigate 
this curious outcome, the IC contribution to the final roll, pitch, and yaw rms errors 
was broken down into the contributions from each of the six different groups of error 
states, as described in Appendix A. This breakdown, given in Table 5.5, shows that the 
final attitude uncertainties, after 10000 seconds, are insensitive to the initial attitude, 
attitude rate, torque, and line bias uncertainties and are only minutely sensitive to the 
initial lever arm conditions. The bulk of the IC contribution is due to multipath, as a 
result of reinitialization when a satellite dropout or switch occurs. Recall that in the 
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state transition matrix, the 4x4 block matrix, corresponding to the channel with the 
switch, is set to identity and all associated off-diagonal elements are zeroed. The last 
switch occurred only 30 seconds from the end of the run, and naturally, the effects of 
that reinitialization are still felt: the process noise (PN) relative contribution decreases, 
while the IC relative contribution increases. The effect on the yaw is larger than on the 
roll and pitch, consistent with the larger dips in yaw uncertainty of Figure 5.6. It is 
unclear why the yaw is apparently more correlated with the multipath than are roll and 
pitch. 


Table 5.5 Initial Condition Contribution Breakdown: Percentages 


Initial Condition 
Contribution 

Per Cent of Total 
Variance 

Roll 

Pitch 

Yaw 

Attitude 
Attitude Rate 
Torque 
Line Bias 

0.00 

0.00 

0.00 

Lever Arm 

0.04 

0.04 

0.05 

Multipath 

1.93 

1.90 

15.04 

Sub-Totals 

1.97 

1.94 

15.09 


Despite the reinitialization, which occurs typically on one channel only, the final 
uncertainties are seen to be dominated by the process noises from the bias-type states - 
the lever arm, multipath, and line bias except for the yaw component. The yaw’s 
insensitivity to line bias may result from the distribution of the satellites in azimuth. A 
rotation could be interpreted as a change in line bias, or vice versa. Concerning yaw, a 
positive change in line bias observed from one satellite would be counteracted by a 
negative change in line bias observed from a satellite 180 deg in azimuth from the 
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first. All observable satellites, however, are in one hemisphere, causing only one- 
signed changes in line bias for roll and pitch rotations. 

Figures 5.16, 5.17, and 5.18 show the IC, PN, and measurement noise contribu- 
tions to the roll, pitch, and yaw variances for the nominal run. Each curve is the per- 
centage of the total variance of that component. The yaw component is easily 
discerned from the roll and pitch, which have very similar profiles. Note the evidence 
of the reinitialization in the multipath contributions. Also note that the line bias IC and 
PN contributions to the yaw component are less than 5% for all time. 
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Figure 5.16 Initial Condition Contributions to Attitude Variances 
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Figure 5.17 Process Noise Contributions to Attitude Variances 
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Figure 5.18 Measurement Noise Contribution to Attitude Variances 
5.3 Off-Nominal Run Results 

A better understanding of system performance requires examination of its sensitiv- 
ity to various filter parameters. This section presents and discusses the results of a few 
of these off-nominal cases, summarized in Table 5.6, for the 49 state optimal filter. The 
“SS Error” for each Markov process is the steady-state, maximum uncertainty used to 
characterize it. 


Table 5.6 Off-Nominal Conditions 


Parameter 

Conditions Considered 

Orbital Inclination 

0, 28.5, 51.6, 90 deg 

Satellite Selection Criterion 

GDOP, ADOP 

Multipath SS Error 

0, 2, 5, 10 mm 

Multipath Correlation 

0, 50, 100% 

Lever Arm SS Error 

0, 2, 5, 10 mm 

Line Bias SS Error 

0, 2, 5, 10 mm 

Torque SS Error 

0, 5x1 O' 7 , 10' 5 , 10' 2 Nm 

Phase Measurement Error 

1, 3, 10, 50 deg 

Attitude Error Rate Process 
Noise PSD 

0, 10' 7 , 10' 2 (deg/s 2 ) 2 /Hz 

Antenna Half-Cone Angle 

70, 85, 95 deg 
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5.3.1 Orbital Inclination 


As discussed in Section 5.2.1, the satellite geometry has a large impact on the atti- 
tude errors at any one time. The orbital inclination is therefore likely to affect the sys- 
tem performance. One might suspect that fewer satellites would be in view over the 
poles. Indeed, fewer satellites are seen at higher elevations, but since the orbital planes 
of the GPS constellation are equally dispersed in longitude, many satellites are visible 
at low elevations with good azimuth distribution when over a pole. A comparison of 
Figure 5.1 and Figure 5.19 actually indicates that GPS satellite visibility is somewhat 
worse for the equatorial orbit than for the polar orbit. Four satellites in view is a more 
common occurrence in the equatorial orbit, occasionally dipping down to three, and 
only once are more than six satellites in view. 



Figure 5.19 Simulated GPS Satellite Visibility in Equatorial Orbit 

Naturally, this degradation in visibility takes its toll on the GDOP and ADOP, as 
well as the attitude uncertainty, though the effect is not a drastic one, as illustrated in 
Table 5.7. Table 5.8 shows the rms of the attitude error magnitude over the last 6000 s 
of each 10000 s run for different inclinations. The nominal results are shaded. There 
seems to be no significant effect of the orbital inclination on the attitude uncertainty 
based on these results. 
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Table 5.7 Attitude Errors Over Final 6000 s in Equatorial Orbit 


Component 

Min RMS 
Error 
(deg) 

Max RMS 
Error 
(deg) 

Average 
RMS Error 
(deg) 

Roll 

0.466 

0.556 

0.519 

Pitch 

0.468 

0.554 

0.515 

Yaw 

0.287 

0.425 

0.344 


Table 5.8 Average Attitude RMS Error Magnitudes vs Orbital Inclination 


Inclination 

(deg) 

Average 
RMS Error 
Magnitude 
(deg) 

0.0 

0.808 

28.5 

0.785 

51.6 

0.773 

90.0 

0.769 


5.3.2 DOP Type 

The figure of merit employed in rating the satellite geometry will obviously impact 
the attitude solution. The issue here is the effect of using GDOP, rather than the more 
appropriate ADOP, on performance. For the nominal case in which an 85 deg half- 
cone angle was used, the difference in errors between the two selection criteria was 
minimal since the number of visible satellites rarely exceeded the maximum number 
of six which can be tracked at one time. (Refer to Figure 5.1 for satellite visibility.) 
The largest difference in ADOP, at 1350 s, corresponded to a 0.1 deg, or 40%, differ- 
ence in attitude error magnitude. To better see the effects of the selection criteria, the 
antenna visibility cone is now opened up with a 95 deg half-angle to provide a greater 
number of visible satellites than in the nominal case. Figure 5.20 shows the average 



















number of satellites in view for this case to be approximately seven, as opposed to five 
in the nominal run. 



Figure 5.20 Simulated GPS Satellite Visibility: 95° Half-Cone Angle 


In Figure 5.21, we now see more of a divergence between the best ADOP and the 
ADOP resulting from the satellites yielding the best GDOP. To translate this ADOP 
improvement into an attitude error improvement, define the difference in error magni- 
tude to be: 


where 



*sei = ADOP sel o v and c best = ADOP bes[ o v 
and the subscripts “sel” and “best” refer to the selected satellites based on the smallest 
GDOP, and the satellites giving the best ADOP, respectively. The per cent difference, 
plotted in Figure 5.22, is defined as: 

Per Cent Diff = —100% 

°sel 

We see that, with this definition of improvement, using the ADOP criterion over 
GDOP results in a maximum decrease in rms error magnitude of 55% for this case. 
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The average difference over the 10000 s shown constitutes a sizable 23% improve- 


ment. 




Figure 5.22 Per Cent Difference in Attitude Error Magnitude: 

GDOP vs ADOP 


The DOP criterion is only an indicator of the effects of instantaneous satellite 
geometry on the states of primary interest, assuming no previous knowledge. It is not a 
predictor of system performance, as implied by Table 5.9, because other states are 
involved and the covariance prior to the measurement is not diagonal and infinite. 
Though the states of primary interest in this application are the roll, pitch, and yaw. 
Table 5.9 actually shows the performance in attitude to be slightly better when using 
the GDOP figure-of-merit. It is possible to choose the set of satellites which minimizes 
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the attitude uncertainty given the existence of the other 46 states and the finite covari- 
ance, but this would defeat the purpose of using the relatively easily calculated figure- 
of-merit. Regardless, the selection criteria will become less of an issue as hardware 
technology improves, allowing more satellites to be tracked at one time. 


Table 5.9 Attitude Performance vs Selection Criterion 


Selection 

Criterion 

Average RMS Errors (deg) 

RSS 

(deg) 

Roll 

Pitch 

Yaw 

GDOP 

0.483 

0.483 

0.289 

0.742 

ADOP 

0.487 

0.485 

0.300 

0.750 


5.3.3 Multipath 

The error budget in Section 5.2.6 indicates that the multipath states have a large 
impact on the attitude error statistics for the nominal run. The attitude error sensitivity 
to these states is examined further here by varying separately the maximum uncer- 
tainty and the correlation between antennae about the nominal values, <5^ = 5 mm and 
= 0.5. 

Table 5.10 shows the maximum roll, pitch, and yaw rms errors and the average 
rms errors over the last 6000 s of the run for = 0, 2, 5, and 10 mm. Not surprisingly, 
the performance substantially degrades with increasing multipath noise. Plots of the 
attitude errors (not provided here) show that the magnitude of the jumps at satellite 
switches also increase with multipath noise and are more pronounced in the yaw com- 
ponent. 

Table 5.11 below shows the attitude error statistics for 0%, 50%, and 100% corre- 
lation of the multipath states between the four antennae to each satellite. Clearly, the 
attitude knowledge becomes better with multipath correlation. Note that the rms errors 
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Table 5.10 Attitude RMS Errors vs Multipath Level 


Multipath 
Level, G? 
(mm) 

Max RMS Errors (deg) 

Average RMS Errors (deg) 

Roll 

Pitch 

Yaw 

Roll 

Pitch 

Yaw 

0 

0.325 

0.323 

0.204 

0.293 

0.294 

0.194 

2 

0.399 

0.397 

0.257 

0.369 

0.367 

0.221 


ggggg 


IS 


0.494 


10 

0.803 

0.787 

0.719 

0.721 

0.722 

0.523 


fc le 100% correlated case are the same as those for = 0, shown in Table 5.10. In 
the 100% correlated case, the filter is effectively only estimating six multipath states, 


one per satellite, as opposed to 24 states for the uncorrelated or partially correlated 


ca _ . The multipath states for each antenna to one satellite are indistinguishable and 

£ 

dr j out in the phase difference, and their uncertainties consequently remain at the 

m t mum 5 mm value, 
w 


5 

w 

C 

Q 

u 

Q 
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Table 5.11 Attitude RMS Errors vs Multipath Correlation 


Multipath 

Correlation 

(%) 

Max RMS Errors (deg) 

Average RMS Errors (deg) 

Roll 

Pitch 

Yaw 

Roll 

Pitch 

Yaw 

0.0 

0.637 

0.632 

0.539 

0.587 

0.587 

0.399 

50.0 

0.536 

0.533 


mm 



100.0 

0.325 

0.323 

0.204 

0.293 

0.294 

0.194 


5.3.4 Lever Arm 

The error budget for the nominal run reveals a strong sensitivity of the attitude 
uncertainty to the lever arm noise level. This sensitivity is confirmed in Table 5.12, 
which provides the attitude uncertainties for variations in the this parameter. The same 
trends are seen: yaw is consistently more observable than the roll or pitch, and all three 
become less observable with poorer knowledge of relative antenna position. 
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Table 5.12 Attitude RMS Errors vs Lever Arm Noise Level 


Lever Arm 
Noise Level, 
<?] (mm) 

Max RMS Errors (deg) 

Average RMS Errors (deg) 

Roll 

Pitch 

Yaw 

Roll 

Pitch 

Yaw 

0 

0.468 

0.460 

0.367 

0.422 

0.418 

0.256 


mm 





| 

5 

0.790 

0.799 

0.600 

0.764 

0.771 

0.529 

10 

1.308 

1.348 

0.998 

1.289 

1.325 

0.951 


5.3.5 Line Bias 

The results of varying the maximum line bias uncertainty, provided in Table 5.13, 
concur with the error budget tabulations: the yaw rms error is fairly insensitive to line 
bias errors. There is only a 7% increase in yaw variance from cp = 0 to Gp = 10 mm, 
compared to a 300% increase in roll and pitch variances. Why the line bias error 
affects roll and pitch more than yaw is not entirely clear, but it may have to do with the 
satellite distribution in azimuth as mentioned in the nominal error budget discussion. 


Table 5.13 Attitude RMS Errors vs Line Bias Noise Level 


Line Bias 
Noise Level, 
Op (mm) 

Max RMS Errors (deg) 

Average RMS Errors (deg) 

Roll 

Pitch 

Yaw 

Roll 

Pitch 

Yaw 

0 

0.460 

0.441 

0.408 

0.399 

0.399 

0.315 


0.536 

Ijlgl 

EBB 



0.316 

5 

0.742 

0.729 

0.430 

0.660 

0.652 

0.321 

10 

0.991 

0.956 

0.450 

0.826 

0.811 

0.326 


5.3.6 Torque 

The attitude uncertainties about all three axes are affected very little by the torque 
noise level, as evidenced by Table 5.14. The nominal value used for the maximum 
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torque uncertainty is 5% of the gravity gradient torque acting on the spacecraft for 
pitch and roll angles of 5 deg. The greatest value, 0.01 Nm, is four orders of magnitude 
larger than the nominal, and yet the attitude variances increase by less than 9%. This 
noise in the torque does, however, affect the attitude error rate statistics. These results 
are reasonable since the torque, or acceleration, is integrated to give the velocity, and 
the phase difference is insensitive to both - they must be inferred from the attitude 
through correlations. 


Table 5.14 Attitude RMS Errors vs Torque Noise Level 


Torque Noise 
Level, a T 
(Nm) 

Max RMS Errors (deg) 

Average RMS Errors (deg) 

Roll 

Pitch 

Yaw 

Roll 

Pitch 

Yaw 

0 

0.536 

0.533 

0.414 

0.496 

0.495 

0.317 


0.5362 

si 

SIllS 

mmm 

MMM 

1|§§§ 


10' 5 

0.541 

0.535 

0.418 

0.500 

0.496 

0.319 

10’ 2 

0.550 

0.543 

0.437 

0.509 

0.505 

0.332 


5.3.7 Phase Measurement Noise 

The effect of the accuracy of the total phase measurement on the attitude perfor- 
mance is also investigated. The results for 1, 3, 10, and 50 deg uncertainties, corre- 
sponding to 0.5, 1.6, 5.3, and 26.4 mm, respectively, are given in Table 5.15. The 
attitude accuracy does not appear to degrade as much as one might expect from a 
seemingly large measurement uncertainty, as indicated in the 50 deg case. Simply 
based on an average ADOP of 2.5, the RSS error from one set of measurements, 
between 12 and 18 phase differences given the satellite visibility shown in Figure 5.1, 
would be in the neighborhood of 3.8 deg, but fortunately, a set of measurements is 
taken once every 6 s. The poor measurement accuracy is counteracted by the sheer 
number of measurements, since probability theory asserts that the uncertainty after n 
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samples will be reduced by a factor of Jn (with no process noise of course). This 
insensitivity is also reflected in the nominal error budget and Figure 5.18. 


Table 5.15 Attitude RMS Errors vs Phase Measurement Noise 


Phase Meas 
RMS Error, 
(deg) 

Max RMS Errors (deg) 

Average RMS Errors (deg) 

Roll 

Pitch 

Yaw 

Roll 

Pitch 

Yaw 

1 

0.529 

0.525 

0.407 

0.490 

0.488 

0.310 

3 





0.495 

0.317 

10 

0.570 

0.572 

0.453 

0.527 

0.528 

0.352 

50 

0.825 

0.855 

0.773 

0.765 

0.782 

0.617 


5.3.8 Attitude Rate Noise 

As mentioned earlier, the filter will weight the measurements more heavily with 
larger process noise. Essentially, the filter does not have as much “memory,” and its 
knowledge of the state decreases to a greater extent between measurements. It there- 
fore can not rely on previous measurements to help lower the variances. When the pro- 
cess noise is large, the best the filter can do is bring the variance down to a level 
dictated by the measurement accuracy and the instantaneous geometry, as indicated by 
the dilution of precision. White noise was added to the attitude rate error state to 
account for unmodelled torque noise. Table 5.16 shows, not surprisingly, the roll, 
pitch, and yaw uncertainties to increase with the white noise level. An interesting 
result does appear, however, in the relative level of uncertainty between the roll and 
pitch: with no noise, roll is more observable than pitch, and adding noise obscures the 
difference. Intuitively, the pitch should be better estimated since any rotation about the 
roll axis is unobservable to the second baseline, which lies along that axis. Recall from 
Section 3.2.2, however, that the roll and yaw are coupled in the equations of motion. 
The roll can therefore be inferred via this coupling. The additive process noise 
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obscures the coupling, effectively reducing the correlation between roll and yaw, and 
the roll thus becomes less observable. 


Table 5.16 Attitude RMS Errors vs Attitude Rate Noise Level 


Q e (deg/s) 2 /s 

Max RMS Errors (deg) 

Average RMS Errors (deg) 

Roll 

Pitch 

Yaw 

Roll 

Pitch 

Yaw 

0 

0.427 

0.491 

0.376 

0.389 

0.445 

0.304 


§gl§j 


lllll 

0.496 



10' 2 

0.555 

0.548 

0.441 

0.513 

0.509 

0.334 


5.3.9 Antenna Half-Cone Angle 

Recall that satellites at lower elevations are more useful in estimating yaw, while 
satellites overhead provide better roll and pitch accuracy. This assertion is supported 
by Table 5.17, which gives the attitude uncertainties for three antenna half-cone 
angles. Since, for RADCAL, a satellite must be in view of all four antennae, restricting 
the visibility cone results in the loss of measurements to satellites at lower elevations. 
Note that the roll and pitch accuracies do degrade with smaller cone angle, and the 
yaw is still more observable. The yaw accuracy, however, suffers more than the other 
two components with tighter viewing restrictions - the average roll and pitch variances 
are onl\ 2 7< ~r greater in the 70 deg case than in the 95 deg case, whereas the average 
yaw variance is 185% greater. 


Table 5.17 Attitude RMS Errors vs Antenna Half-Cone Angle 


Ant Half- 
Cone Angle 
(deg) 

Max RMS Errors (deg) 

Average RMS Errors (deg) 

Roll 

Pitch 

Yaw 

Roll 

Pitch 

Yaw 

70 

0.783 

0.711 

1.078 

0.541 

0.544 

0.488 

: 

EUS 

mmammaesimm 

|jgg| 


wmm 

. 

wmmmmasmm 


95 

0.534 

0.531 

0.370 

0.483 

0.483 

0.289 
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5.4 Optimal-Suboptimal Filters 


It is desirable to reduce the order of the filter from 49 states to something more 
manageable for an onboard navigator. Considered here are two optimal-suboptimal fil- 
ters (using consider state analysis): a 9 state filter which estimates attitude, attitude 
rate, and torque, and a 13 state filter which includes the 4 line biases. (If the double 
difference measurement is used instead of the single difference, the line bias states 
could potentially be removed from the filter since they are specific to the antennae and 
therefore cancel. The performance in using double differences, however, would be 
worse, since this is equivalent to assuming infinite line bias variance in the filter using 
single differences; the filter is not taking advantage of the limits on line bias uncer- 
tainty [12].) As previously discussed, the lever arm and multipath states (altogether, 36 
= 12 lever arm + 24 multipath) appear to have a significant impact on system perfor- 
mance errors. The effect of their omission is examined by “considering” them as 
described in Section 4.1.3. The results presented here are the attitude errors from these 
so-called “optimal-suboptimal” filters. 

Figure 5.23 shows the time histories of the attitude rms errors and the magnitudes 
(the RSS of the roll, pitch and yaw errors) for the two optimal-suboptimal filters as 
well as the full optimal filter. The same initial conditions and process noise parameters 
are used in all three cases. Not surprisingly, the 9 state filter’s performance is the 
worst. The 13 state filter is seen to perform somewhat better, in general. 

The results are summarized in Table 5.18, which shows the average attitude rms 
errors and their RSS’s over the last 6000 s of the 10000 s run. The results of the nomi- 
nal error budget (Section 5.2.6) might lead one to believe that a filter which does not 
estimate the lever arm, line bias, and multipath states would perform quite poorly 
compared to the optimal filter. The table shows, however, that the attitude’s average 
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RSS Attitude Error (deg) Yaw Error (deg) Pitch Error (deg) Roll Error (deg) 




Figure 5.23 Suboptimal Filter Performance 
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mean square errors of the 9 state filter are less than 26% larger than the 49 state “truth” 
filter. The average attitude variances of the 13-state filter are less than 18% larger than 
those of the optimal filter. This seems to imply that the optimal filter does not do as 
good a job at estimating the lever arm, line bias and multipath states as might be 
thought. 


Table 5.18 Optimal-Suboptimal Filter Performance 


Filter 

Average RMS Errors (deg) 

RSS 

Roll 

Pitch 

Yaw 

(deg) 

9 State O-S 

0.552 

0.555 

0.346 

0.856 

13 State O-S 

0.528 

0.534 

0.344 

0.826 

Optimal 

0.496 

0.495 

0.317 

- 

0-769 


The error budgets for the 9 state and 13 state suboptimal filters are provided in 
Table 5.19 and Table 5.20, respectively. A comparison with the error budget of the 
nominal run reveals that the IC contribution, due essentially to the multipath IC, is 
larger for the optimal-suboptimal filter. Note also that, for the 9 state filter, the line bias 
process noise contributes more to the attitude uncertainty, in an absolute as well as rel- 
ative sense, than in the optimal case (e.g., the line bias contribution to roll error: 
35.18%, 0.298 deg vs 23.16%, 0.221 deg). In adding the line bias states, their contri- 
butions are seen to decrease in Table 5.20. 

One might expect the contributions from the states which are not estimated to be 
larger, in an absolute sense, than in the optimal filter. This is indeed the case in the line 
bias PN contribution for the 9 state filter. Note, however, that the lever arm PN contri- 
bution is only slightly larger, and the multipath PN contribution is actually somewhat 
smaller! The contributions still add to give larger final variances. A crucial point is the 
fact that the optimal filter minimizes the trace of the entire 49x49 covariance matrix. 
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Table 5.19 Error Budget for 9 State Optimal-Suboptimal Filter 



Per Cent of Total 
Variance 

Final RMS Errors 
(deg) 

Contribution 

Roll 

Pitch 

Yaw 

Roll 

Pitch 

Yaw 

1 ) Initial Conditions 

4.08 

3.92 

27.15 

0.102 

0.100 

0.172 

2) Process Noise: 



Attitude Rate 

1.11 

1.34 

4.88 

0.053 

0.058 

0.073 

Torque 

0.00 

0.00 

0.12 

0.002 

0.003 

0.011 

Lever Arm 

26.60 

26.56 

33.86 

0.260 

0.260 

0.192 

Line Bias 

35.18 

35.98 

1.54 

0.298 

0.303 

0.041 

Multipath 

32.08 

31.32 

30.85 

0.285 

0.283 

0.183 

3) Measurement 
Noise 

0.95 

0.87 

1.58 

0.049 

0.047 

0.041 

Totals 

100.0 

100.0 

99.98 

0.503 

0.505 

0.329 


Table 5.20 Error Budget for 13 State Optimal-Suboptimal Filter 



Per Cent of Total 
Variance 

Final RMS Errors 
(deg) 

Contribution 

Roll 

Pitch 

Yaw 

Roll 

Pitch 

Yaw 

1 ) Initial Conditions 

6.54 

6.84 

28.15 

0.126 

0.130 

0.174 

2) Process Noise: 



Attitude Rate 

1.30 

1.21 

3.96 

0.056 

0.055 

0.065 

Torque 

0.00 

0.00 

0.09 

0.003 

0.002 

0.010 

Lever Arm 

28.01 

27.80 

34.11 

0.260 

0.261 

0.192 

Line Bias 

17.92 

17.87 

0.69 

0.208 

0.210 

0.027 

Multipath 

42.45 

42.86 

31.20 

0.320 

0.325 

0.183 

3) Measurement 
Noise 

1.56 

1.51 

1.65 

0.061 

0.061 

0.042 

Totals 

97.79 

98.09 

99.85 

0.491 

0.496 

0.328 


























































































































not the roll, pitch, and yaw variances. Therefore, it is possible for a suboptimal filter to 
estimate some individual states better than the optimal filter, but overall, the perfor- 
mance will be worse. 

Based on this cursory analysis, the 13 state filter appears to be a good choice for an 
onboard system. It provides substantially smaller computer burden with little signifi- 
cant reduction in attitude accuracy. Of course, the performance shown here is the ideal 
overall attitude performance of the reduced-order filters, because the filters considered 
the omitted states (but did not estimate them). The actual reduced-order filter would 
require “tuning” to approach this “best” performance. Filter tuning can involve such 
methods as underweighting and adjusting the process noise. Better attitude accuracy, 
with little additional computer burden, may be also be achieved by lumping the omit- 
ted states into one new state. 
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Chapter 6 


Conclusions 

The goals of this thesis were to determine the attitude accuracy achievable in pro- 
cessing GPS phase difference measurements with a Kalman filter and to understand 
the effects of the major error sources on that accuracy. The performance of a GPS 
interferometric attitude determination system was investigated via a linear covariance 
analysis. Error sources were identified and the significant contributors were modelled. 
The optimal system’s sensitivity to these error sources was examined through the error 
budget and by varying system parameters. The performance of two optimal-subopti- 
mal filters was also investigated as a first step toward designing an onboard filter. A 
summary of the results and suggestions for future work are provided here. 

6.1 Summary of Results 

For all runs investigated, the yaw component of the attitude consistently exhibited 
markedly different behavior from roll and pitch. In fact, the yaw was repeatedly better 
estimated than roll and pitch, due to the availability of lower-elevation satellites which 
provide larger yaw-sensitivity of the phase difference measurement. This sensitivity 
was observed in varying the antenna field-of-view: the yaw suffered more than the roll 
and pitch with decreasing cone angle. Under nominal conditions, the average rms 
errors for roll, pitch, and yaw were 0.496, 0.495, and 0.317 deg, respectively. In the 
presence of attitude rate noise, roll and pitch were observable to the same degree, 
while in its absence, roll was more observable because no noise obscured its coupling 
with the better estimated yaw component. Yaw was found to be fairly insensitive to 
uncertainties in the electrical path lengths to the antennae, compared to the other two 
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components. Lastly, yaw showed more correlation with the multipath states, as evi- 
denced by the more pronounced changes (compared to roll and yaw) in its error with 
satellite switches. 

After close inspection, it was found that the results of the simulation are consistent 
with the multipath model, which assumed that a new satellite is uncorrelated with the 
spacecraft’s state. In truth, the multipath, though complex in nature and difficult to 
characterize, behaves in a predictable way and depends on the local reflective environ- 
ment and the direction of signal incidence. In other words, the multipath from a newly 
introduced satellite should, in fact, be correlated with the attitude. 

Based on the results of the error budget and the off-nominal runs, it was seen that 
good knowledge of relative antenna locations, path delays, and the multipath environ- 
ment is required for high attitude accuracy. The error budget of the nominal run 
showed that roughly 70% of the roll, pitch, and yaw variances was due to noise in the 
lever arm and multipath states. (The lever arm flexures, however, are probably more 
easily controlled, i.e., more easily calibrated out, than multipath.) The noise in the line 
biases were seen to contribute about 23% to the roll and pitch variances, while contrib- 
uting a meager 1% to the yaw variance. These results were confirmed in varying the 
error levels of these states. Performance improved significantly in all three compo- 
nents when the lever arm flexures and the multipath levels were reduced. Better accu- 
racy was also achieved with greater multipath correlation between antennae. And 
though yaw accuracy remained largely unchanged, roll and pitch were better estimated 
with smaller line bias uncertainty. 

Several factors which were thought to play a significant role in attitude perfor- 
mance but in fact did not were the torque error level, the measurement uncertainty, 
orbital inclination, and satellite selection algorithm. High accuracy was maintained 
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with large torque uncertainties. When the torque error level, a T , was increased by four 
orders of magnitude from its nominal value, the attitude variances increased by less 
than 9%. The performance did not considerably degrade with seemingly poor mea- 
surements. Yaw was seen to experience the largest deterioration; even so, its average 
variance increased by only 23% with a change in total phase measurement uncertainty 
from the nominal 3 deg to 10 deg. These insensitivities were also reflected in the error 
budget. The accuracy actually improved slightly with larger orbital inclinations. 
Finally, virtually no difference in performance was observed between selection crite- 
ria. 

In summary, the results of the nominal and off-nominal cases indicate that attitude 
accuracy is limited by: 

• Lever arm errors 

• Line bias errors 

• Multipath errors and correlation 

• Antenna field-of-view 

w hile it is fairly insensitive to: 

• Satellite selection criterion (GDOP vs ADOP) 

• Orbital inclination 

• Torque errors 

• Measurement uncertainty 

The statistical performance of two suboptimal filters, using consider states, was 
also investigated. The first estimated 9 states, consisting of the attitude errors, their 
rates, and the torque error, and the second included the 4 line bias error states. In con- 
sidering the apparently major contributors to attitude errors, the suboptimal filtering 
results did not show as large a deterioration in performance as the error budget and off- 
nominal runs implied. The average attitude variances in the 9 and 13 state filters were 
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only 26% and 18% larger than the 49 state optimal filter. This result implies that the 
optimal filter is not estimating the lever arm, line bias and multipath error states as 
well as might have been originally thought. Consistent with the nominal error budget, 
however, were the relative improvements in the attitude variances between the 9 and 
13 state filters. Based on the optimal error budget, the yaw component is insensitive to 
line bias errors, and indeed, the yaw accuracy improved very little in adding those 4 
states while the roll and pitch variances decreased by 10%. It is important to remember 
in comparing results of suboptimal and optimal filters, however, that the optimal filter 
minimizes the trace of the entire covariance matrix, not the attitude variances. It is 
therefore possible, though unlikely, that a suboptimal filter could be computed to give 
better attitude performance than the optimal filter. But based on this cursory analysis, 
the 13 state optimal-suboptimal filter represents a good compromise between accuracy 
and computational burden. 

6.2 Suggestions for Future Work 

A logical first step in continuing the investigation of GPS interferometry applied to 
the attitude determination problem considered here would be to implement models 
which better characterize the error sources. As previously mentioned, the lever arm 
and line bias error models assumed in this study attempt to mimic the effect of thermal 
variations. An improvement in their characterizations would be to take into account 
the correlations which arise from the physical phenomenon assumed to drive these 
errors. The lever arm flexures are likely to behave similarly (as opposed to their 
assumed independence) when undergoing heating, for example. Another measure 
which would yield more realistic performance is to include inertia matrix errors in the 
state. This is because body axes and principal axes are always somewhat different. The 
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multipath model could be improved by accounting for its correlation with the attitude 
when a new satellite is introduced. One possible amendment is to maintain the correla- 
tions with the multipath states associated with each satellite. This could be achieved 
by including the multipath states associated with all 24 operational GPS satellites (thus 
increasing the state size from 49 to 121 !). Finally, though the attitude performance was 
found to be largely insensitive to torque errors, the effects of solar radiation, atmo- 
spheric drag, and other sources of torques could be modelled more accurately. Also, a 
more in-depth analysis would include an investigation of the performance in off- 
LVLH nominal attitudes. 

The next step in investigating system performance would be to implement a Monte 
Carlo simulation. The Monte Carlo simulation involves determining the statistical per- 
formance empirically through multiple runs and actual state estimation. This would be 
useful in verifying the validity of the covariance analysis results. Nonlinear effects 
such as integer ambiguities and cycle slips could be included to enhance the fidelity of 
the model. 

Another area of further research is in suboptimal filtering. The performances of 
two optimal-suboptimal filters were summarily examined in this study. A more in- 
depth examination would involve the analyses of filters with other combinations of 
states and their tunings. Various methods of filter tuning could be explored, such as 
adjusting the process noise levels and underweighting the measurements. Another tun- 
ing technique is to approximate the observed optimal-suboptimal gains with analytic 
functions of time, typically exponentials and piecewise constants [10]. The result 
would be a filter suitable for onboard navigation. 

Finally, actual flight data can be filtered, since that is the ultimate goal for the 
aforementioned exercises. As part of this, the whole area of initial attitude determina- 
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tion for an onboard system with multipath errors needs to be addressed. RADCAL 
provides one platform for investigation, and in the near future, other opportunities will 
arise in the increasing number of applications of GPS-based attitude determination 
systems. 
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Appendix A 


Error Budget Implementation Issues 

The major steps in the tabulation of the error budget are shown in the simulation 
flow chart in Figure 4.4. The measurement update and propagation steps for the error 
budget differ, however, from those of the total covariance. These steps are described in 
greater detail here. 


A.l Measurement Update 

At the measurement time, the best satellites are first selected. Just as the multipath 
elements of the total covariance matrix must be reinitialized when a satellite has 
dropped out of view, so must the contributions if the error budget is desired. To effect 
this reinitialization in the contributions, the state transition matrix, <I>, is reset by insert- 
ing a (4x4) identity matrix in the block corresponding to the lost satellite and zeroing 
the associated off-diagonal terms. The process noise constituent is further decomposed 
into the contributions from the noise on the attitude rate, torque, lever arm, line bias 
and multipath states: 



+ P e T + P G, + P G« 



The process and measurement noise contributions are reinitialized by simply zeroing 


out the entire block row and block column corresponding to the appropriate channel. 


With the satellites selected and the multipath states in each component updated, 
the measurement sensitivity and thus the weighting, optimal or otherwise, is calcu- 
lated. The measurement is then incorporated into the covariance components via equa- 
tions (4.22). 
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A.2 Propagation 

In propagating the total covariance matrix, recall that the matrix Ricatti equation is 
integrated using a fourth-order Runge-Kutta scheme. For the error budget, however, 
we opt for the discrete form of the propagation which employs the state transition 
matrix as opposed to the dynamics matrix, and the process noise covariance instead of 
its power spectral density. The incremental state transition matrix, cp , can be calcu- 
lated by simply integrating the dynamics matrix, F, over the time step. This is done 
using the Runge-Kutta integrator that is already in place. The process noise covari- 
ance, N, is not as easily attained for the error budget calculations given in equations 
(4.20). N is decomposed into its components for each state to propagate each of the 
aforementioned process noise constituents. An approximation to it is found by expand- 
ing Equation (4.19), rewritten here: 

N = <p(t k , X)Q<p T (t k ,T) dx 

J h-, 

in a truncated Taylor series. The derivative of N is found by applying Leibnitz’s rule to 
the above equation [ 20 ]: 


N = y ^-[cp(^t)Qcp r (^,t)]rfx+ [cp(/ t ,f t )Q<p r (f*,f*)]^- (A.i) 


h-\ ul k 


+ [9 (**> t k-\) (**> h- 1 )] ^ t 

dt k _ . 

for infinitesimal A t=t k - t k .j. Noting that - 3 — = 0 and that <p (t k , t.) = I , and tak- 
a'* 

ing Q to be constant over Ar, the time derivative of N reduces to: 

N = Q + F ^ <p(f t , T)Q(p\t k , x) dx+ <p(f 4 , x)Q(p\t k ,x) dxF T (A.2) 

Jl k-i J ' k -i 

= Q + FN + NF 7 (A.3) 

Assuming F to be constant over A t, the second and third derivatives of N are obtained: 
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(A.4) 


N = FN + N F r 


N = FN +NF r (A.5) 

These expressions for the derivatives are inserted into the Taylor series, which is: 

N = N 0 + N 0 Af + ^N 0 At 2 + j,N 0 Ar 3 + ... (A.6) 

where N 0 is the nominal noise covariance, i.e., the covariance at time t k .j\ thus, N 0 = 0 
and N o = Q. Since the process noise contributions to the decomposed total covariance 
matrix are each handled separately, the first derivative of the 49x49 noise covariance 
of process, x, is: 


N 0 = 

X 


0 0 

Q 


X 

0 0 


ii 49x49 


where Q x is the appropriately dimensioned noise PSD associated with process x. The 
higher derivatives are then computed using the assumed constant dynamics matrix, F 
in equations (A.4) and (A.5) and inserted into the series in (A.6), truncated to four 
terms. (Four terms were deemed sufficient, since over 10000 s the maximum differ- 
ence in the trace of P between the two methods of propagation - integrating the Ricatti 
equation and using the state transition formulation - was less than 0.05%.) The state 
transition formulation is then used in propagating each noise component: 


P q x = ♦P flj * 7 +N x 

A.3 Final Tabulation 

At the end of the run, the covariance breakdown is now tabulated. The breakdown 
in full is: 


P„ , = <t>p A o 

r Total 


+ P, 


+ P 2 T + P 2. 


+ P n +P n +P 


M 
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Extracted from each piece of the total covariance are the first three diagonal terms - the 
contributions to the attitude variances. The attitude variances in the first term on the 
right arise from all states’ initial conditions. This piece can be decomposed even fur- 
ther to see the contributions from the initial conditions of each group of states at no 
extra cost: 

<pp.<p r = <b [P 0 + P a + P T + P, + P R + P £ ] <D r 

0 °0 00 *0 *0 Po So J 

In performing the error budget, a total of eight 49x49 matrices - the total covari- 
ance, the state transition matrix, the five process noise covariance matrices, and the 
measurement noise contribution - must be propagated and updated with measure- 
ments. Though the number of operations can be reduced somewhat by utilizing the 
fact that the dynamics and process noise matrices are sparse, the error budget is still a 
costly undertaking, but a very illustrative one. 
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Appendix B 


Discussion of RADCAL Flight Data 

Several hours of data from RADCAL’s TANS Quadrex receiver was obtained from 
Penny Saunders at Johnson Space Center. The data is in binary form and is parsed into 
readable ASCII format using a program called TANSPOST, also provided by JSC. In 
addition to furnishing the raw differential phase data, TANSPOST gives RADCAL 
position fixes, velocities, latitudes, longitudes, and altitudes, satellite selection infor- 
mation (PDOP, HDOP, VDOP, TDOP), and receiver health information. The TANS 
receiver computes position solutions from the GPS Coarse/ Acquisition (C/A) code, 
subject to selective availability (SA). Virtually no accuracy is lost in the attitude solu- 
tion due to SA, however, since the lines-of-sight are insensitive to the position errors 
induced by the intentional degradation. 

The raw phase difference measurement files contain phase differences only to sat- 
ellites in view of all four antennae, even though a measurement could have been pro- 
cessed if only two antennae view a particular satellite. The RADCAL system is 
nevertheless still overdetermined given the 18 measurements to the three baselines. 
The data files obtained for this investigation corresponds to GPS week 708, day 6, or 
August 7, 1993. Figure B.l shows the PRN’s which are tracked for a 10000 s interval, 
and Figure B.2 shows the number of satellites in RADCAL’s best view for the same 
interval. Note that the maximum number is six because there are only six channels per 
antenna, and the receiver requires that a satellite must be in view of all four antennae 
to take phase differences. The figure shows that at least four satellites will be in view 
90% of the time. The 85 deg antenna half-cone angle in the simulation was chosen to 
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Figure B.2 Number of Satellites in RADCAL’s Best View 


The phase differences to a satellite are referred to a master antenna, so denoted in 


the raw data file. The master antenna can switch for each satellite based on the signal- 
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to-noise ratio (SNR). It is therefore useful to designate a constant pseudo-master 
antenna, and recompute the phase differences to refer them to this pseudo-master. 
(Note the antenna indices in Figure 1.2. Antenna 1 is the pseudo-master.) For example, 
if antenna 2 is the master, we can refer the phase differences to antenna 1 by: 

A *f > l — » 3 = ^2 — > 3 ~ A< J > 2 — > 1 

A <t>1^4 = A ^2 — » 4 — A 4^2 — > 1 

where A<J> f - _* • is the phase difference across the baseline from master, i, to slave, j. The 
actual phase differences to PRN 19 for each of the three baselines over a 10000 s inter- 
val of time were calculated and are shown in Figure B.3. The outages and cycle slips 
are painfully apparent. In the actual RADCAL flight data, we see phase differences 
which can be greater in magnitude than the baseline length. This results when the 
receiver loses lock and reacquires, and then differences the erroneous phases. It is pos- 
sible then to have a phase difference of thousands of cycles, for different receivers 
[12]. Fortunately, this case involves only one receiver on a small platform, keeping the 
phase predictions during reacquisition comparable. It is also possible that some pre- 
processing on the phase differences was performed, either by the receiver or in TANS- 
POST. Except for the unusually long interval from 1000 to 2800 s in baseline 2 in 
which no slips occur, these plots are typical of the phase differences for other satel- 
lites. 

It may be desirable to pre-process the phase differences even further by removing 
the cycle slips and outlying points. Of course, an onboard filter would have to deal 
with these situations in real time. The cycle slips are ideally integer jumps, but the 
RADCAL data can exhibit non-integer jumps. To remove the slip and patch the phase 
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Figure B.3 RADCAL Phase Difference Measurements to PRN 19 


difference segments together, the time interval between measurements over which the 
slip occurs must be defined. Depending on the time interval, the dynamics between 
measurements can cause the phase difference to change enough to cause a slip to have 
an ambiguous number of cycles. 

The noise on the phase differences can be characterized by curve fitting, via the 
method of least squares, a linear segment and noting the mean square error. A simple 
second-order polynomial fitted to the linear segment of the phase difference for base- 
line 2 in the above plot between t = 2000 s and t = 2800 s gives a change in phase dif- 
ference of 0.0291 +/- 0.0131 cycles (5.5 +/- 2.5 mm) over the nominal 6 s time 
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interval. If the slip occurs over a short enough interval, the jump may be identified as 
an integer. For longer intervals, the gap may be bridged by using a simple curve fit to 
extrapolate the phase differences on either side to the center of the interval [12]. Hope- 
fully, the jump there is sufficiently close to an integer. For long outages which pre- 
clude this extrapolation technique, the integer must be resolved as in the initialization 
process, but by this time, the attitude should be known well enough to easily determine 
the ambiguity. 

The first attempt made at attitude initialization used a simple, brute-force 
approach. The first minute of data, over which no cycle slips occurred and the dynam- 
ics are essentially constant, was used. The single phase differences were differenced 
again in time; the integer ambiguities (one per satellite per baseline) and constant part 
of the line biases (one per baseline) therefore dropped out. The phase differences were 
predicted based on an assumed initial attitude and were differenced in time, like the 
real data. The best attitude should minimize the following cost function: 

• /<0) = S I[(a*;-a^)-(a>;-a >;)] 2 

i = 2 j = 1 

where i indexes over the epochs, j indexes over the three baselines, and the carat 
denotes the predicted value. The above cost function was computed for the range of 
initial attitudes in which roll and pitch were between -20 and 20 deg and yaw was 
between 0 and 360 deg, at 5 deg increments (9x9x72 = 5832 initial guesses). This 
search did not result in a clear-cut initial attitude, due most likely to the noise and the 
multipath on the actual phase differences. 

A second, more sophisticated, attempt at initialization was made [24], In this tech- 
nique, the integers were varied to arrive at possible solutions. For each baseline, two 
between-satellite double difference were formed, and the constraint of the baseline 
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length was used to deterministically solve for each baseline vector in terms of the dif- 
ferences of the lines-of-sight and their cross-product, to within a sign ambiguity. 
Cycling through the allowable integers for each baseline gives multiple solutions, the 
number of which can be reduced by constraining the relative baseline geometry. 
Unreasonable solutions, such as one indicating the craft to be upside-down, can be 
ruled out as well. This still yielded several solutions. The same procedure was per- 
formed in using a different satellite at the same epoch, but unfortunately, the two sets 
of initial attitudes had no solutions in common. Another alternative would be to use 
Knight’s method [14], which would result in a single most likely solution, though that 
solution may not be correct due to RADCAL’s multipath environment. Some sort of 
batch estimation which accounts for the dynamics, as suggested by Axelrad and Ward 
[2, 28], is most likely required to solve for the initial attitude. 
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